Mechanical characterization and induced crystallization in nanocomposites of thermoplastics and carbon nanotubes

Nanocomposites built from polymers and carbon nanotubes (CNTs) are a promising class of materials. Computer modeling can provide nanoscale views of the polymer–CNT interface, which are much needed to foster the manufacturing and development of such materials. However, setting up periodic nanocomposite models is a challenging task. Here we propose a computational workflow based on Molecular Dynamics simulations. We demonstrate its capabilities and showcase its applications, focusing on two existing nanocomposite materials: polystyrene (PS) with CNT and polyether ether ketone with CNT. The models provide insights into the polymer crystallization inside CNTs. Furthermore, the PS+CNT nanocomposite models are mechanically tested and able to predict an enhancement in Young’s modulus due to the addition of highly dispersed CNTs. We accompany those results with experimental tests and provide a prediction model based on Dynamic Quantized Fracture Mechanics theory. Our study proposes representative simulations of polymer–CNT nanocomposites as promising tools to guide the rational design of this class of materials.


INTRODUCTION
One of the most challenging tasks in materials science is the rational design of nano-composites; which are materials that comprise different components, their dimensions being within the nanoscale range.Among the most promising ones are the polymeric nano-composites (PNCs) 1,2 that combine carbon allotropes such as graphene layers (GLs) and carbon nanotubes (CNTs) with plastic polymers.GLs and CNTs have high strength, while polymers can bear large deformations; and it is expected that by controlling their assembly, we should also be able to fine tune their mechanical properties 3,4 .Although there is a growing demand for PNCs 5 , their improvement faces a number of challenges.On the one hand, there are manufacturing problems, such as CNT aggregation and harsh mixing procedures 6,7 .On the other hand, there is an incomplete understanding about the internal nano-composite structure.Polymeric entanglements and the high surface area of CNTs dominate the molecular interplay, affecting the mechanical properties.Therefore, an atomic level description is central for further development of PNC applications.
Molecular dynamics (MD) simulations are well suited to reveal atomic details about PNCs and their interactions.Initially, MD studies of polymer-CNT systems were focused on the adsorption of individual or highly disperse polymers on CNTs [8][9][10][11] .Studying the mechanical properties of PNCs instead requires to model a compact polymer block with embedded CNT reinforcements, a task that is not straightforward.In this regard, a variety of nanocomposite compact models have been built; for example, polyethylene with CNTs [12][13][14] , epoxy with CNTs 15,16 , and polymethyl methacrylate with CNTs 17 .Those nanocomposite models were built with different protocols that usually start with a disperse CNT distribution, then low-density polymers are filled in the empty space surrounding the CNTs, and finally the components are mixed with a MD annealing cycle.In general, building an atomic model of a PNC should require two steps: (1)  filling a periodic box with a compact and well-equilibrated polymer melt, and (2) inserting the reinforcement molecules within the periodic melt.Regarding the first step, several approaches have been proposed; such as Monte Carlo steps 18,19 , where monomers are linked into a growing polymer until box saturation; Hierarchical Coarse-Grained techniques 20 , where softblobs that represent polymer segments are gradually backmapped into monomers and all-atom models; and MD protocols that involve Annealing Cycles 21,22 , where polymers are arranged into a grid and then collapsed through MD simulations at different temperatures and pressures.While Monte Carlo and Hierarchical Coarse-Grained approaches efficiently reduce the computational cost of shoehorning polymers within a finite volume; those methods are not exempt of caveats: Monte Carlo insertions progressively decrease the accessible space and end up generating voids, while Hierarchical Coarse-Grained requires the availability of accurate coarse models as the blobs get close to the allatom scale.Annealing Cycles rely only on all-atom MD simulations; thus, it circumvents the previous limitations but at high computational cost.Still, a problem arises when Annealing Cycle protocols are applied to polymers of high molecular weight: the longer the polymer, the larger the grid spacing; hence, polymers have considerable empty space around them and can self-fold before contacting each other.Regarding the second step, a timeweighted relaxation algorithm, FADE, has been successfully used to insert bead-polymers into melts and buckyballs into argon boxes 23 .FADE inserts a molecule with the intermolecular forces turned off and then gradually turn them on.However, due to the molecular topology of CNTs, FADE cannot be used to embed CNTs into amorphous phases containing long polymers, as the benzene rings would appear in the middle of covalent bonds, resulting in polymer chains crossing the aromatic planes.Similar efforts have been focused on inserting proteins into lipid bilayers 24 ; e. g., removing lipid molecules and using repulsive forces to create a hole 25 or radially growing a protein perpendicular to the membrane plane 26 .Nonetheless, these approaches were designed for small molecules and cannot be extended for long polymers.For example; lipids tails are aligned and about 1 nm long, thus the removal of lipid molecules does create a local cavity.Conversely, long polymers penetrate the periodic cell in all directions and removing polymers creates a low-density region around the cavity.
Here, we present an Annealing Cycle workflow to build nanocomposites of increasing complexity.The workflow steps are included in a computational tool called "nanocomposite", which works as a plugin in the visualization software VMD 27 .Below, we describe the Annealing Cycle workflow and showcase its applications, focusing on two existing PNC materials: polystyrene (PS) with CNT 28 and polyether ether ketone (PEEK) with CNT 29,30 (Fig. 1).First, we build amorphous polymeric melts.Second, we insert double-walled CNTs (DW-CNTs) into the polymer melts.Our results give insights into the CNT-polymer interface, allowing us to observe polymer crystallization.Third, we mechanically test the atomic models comprised of PS and DW-CNTs and report an enhancement of Young's modulus due to DW-CNT insertion.Moreover, we accompany these MD simulations with experimental tests, and present a scaling rule for the dependence of the mechanical properties with loading velocities.

RESULTS AND DISCUSSION
Building polymer melts To create models that represent amorphous materials, we have modified a previous Annealing Cycle protocol, so-called collapsing-annealing 21 , which was designed for small polymers.Figure 2a shows the updated workflow applied to build a periodic PEEK melt.The procedure is repeated for PEEK and PS.It starts with polymers of 100 monomer length that are randomly oriented in a 5 × 5 × 5 grid (Fig. 2b).In such arrangement, long polymers have considerable empty space around them and self-fold before arriving to the central region.Therefore, the initial conformation is transformed from a disperse grid into a crowded blob, where polymers are close to each other but separated by a 1 nm thick layer of empty space.The crowded arrangement is generated by (1) calculating the empty space surrounding each polymer, (2) moving the polymer within that empty space towards the grid center, (3) rotating the polymer if it were blocked and iteratively repeating (1) and ( 2) until there is no displacements.We called this algorithm geometric collapse to highlight its differences from the force-collapse aggregation: neither energies nor forces evaluation are calculated; only nearest-neighbor searches, translations and rotations.Figure 2c shows the polymeric grid after the geometric collapse.Note the difference in system sizes (scale bar applies to Fig. 2b, c).Details of the geometric-collapse algorithm are included in Supplementary Methods 2.
The blob arrangement is then collapsed by a force-driven MD simulation at 1000 K.A box is defined at the center of the polymer blob and a radial force towards the center is applied to all carbon atoms located outside the box (see "Methods").The close proximity among polymers favors intertwinding whereas the high temperature hinders self-aggregation, generating a well-mixed polymeric ball.Then, the aggregated polymers are inserted into a periodic box (Fig. 2d) and MD equilibrated at 1000 K in NPT conditions to fill all empty spaces with polymer.In the previous collapsing-annealing protocol 21 , the high temperature liquefied the short polymers and the system reached uniform mass distribution.However, due to the high entanglement of long polymers, relaxation does not occur on sufficiently short timescales.To redistribute the mass at the boundaries of the periodic box, we now include a expansion/contraction step (Fig. 2e, f); which performs MD simulations with either an outward or an inward force towards the cell center.Once uniform mass distribution is achieved, we apply an annealing cycle: first 700 K for 5 ns, then cooling steps of −100 K/2 ns, finishing with an equilibration at 300 K until density convergence (Fig. 2 g), which required 130 ns for PS and 210 ns for PEEK.Supplementary Table 1 summarizes these MD simulations.
Our polymer melting started at 1000 K, a range of temperature that is usually applied to randomize of ionic crystal models 31,32 ; however, for a polymer model, it generates stretchings and collisions that require multiple MD minimizations.We noted that a temperature of 700 K is high enough to keep the polymer fluid and decrease the number of minimization restarts.Still, 700 K is a high temperature and its usage has advantages and disadvantages: on the one hand, it allows each polymer to overcome the barriers of a rough energy landscape and obtain a different conformation; on the other hand, it also permits to jump barriers that should not be crossed; namely, chiral carbons can be flipped between R and S enantiomers, and dihedrals angles can change between cis and trans conformations.These problems do not apply to PEEK and PS melts, as there are neither chiral centers nor preference for cis or trans.However, in the case of molecules with precise stereochemistry, such as biomolecules, the structures have to be corrected of stereochemical errors 33 .In Supplementary Fig. 1, we include snapshots of a biological composite built with a variant of this workflow; namely a spider silk model that combines crystalline beta sheets and amorphous regions of protein.
Our criterion for equilibration is constant density 34 , which is a macroscopic characteristic that quickly converges in MD simulations of polymeric melts 35 .The densities are monitored throughout the protocol and reach constant values without significant fluctuations during the final step of annealing at 300 K.For PS melt at 300 K (Fig. 3a), the density increases quickly in the first 15 ns, followed by a slow increase during the next 85 ns.A plateau appears in the last 30 ns, with a value of 1.002 g cm −3 ; that is about 5% difference to the reported density for amorphous atactic PS of 1.047 g cm −3 36 .The PEEK melt follows a similar trend (Fig. 3b), reaching a plateau in the last 50 ns of a 210 ns MD simulation at 300 K. Our simulated PEEK bulk system has a density of 1.193 g cm −3 that differs about 6% from the experimental density of 1.264 g cm −3 for amorphous PEEK with zero crystallinity 37 .These bulk models are used to build the subsequent PNC systems.In the Supplementary Discussion, we include further structural analysis that shows the equilibration of our polymeric melts and nanocomposite models; namely an evaluation of homogeneous density within the systems and quantification of the polymer coiling from radius of gyration calculations and dihedral analysis.

Building PNCs
The basic idea behind the insertion workflow is to open a space with the shape of molecule by pushing away the polymers; then, inserting the reinforcement molecule; and finally annealing the combined structure.Figure 4 shows the workflow for adding a DW-CNT into a periodic PEEK melt.Grid-SMD forces 38 are used to open an empty space.This approach is similar to the Molecular Dynamics Flexible Fitting (MDFF) 39 technique, where biomolecules are force-driven into the volume of a cryo-electron microscopy (cryo-EM) grid.Here, we create a exclusion grid that has the contour of the DW-CNT but inverse the direction of the force, driving the polymers out of the grid region.First, the DW-CNT structure (Fig. 4b) is used to map a grid (Fig. 4c and "Methods").After that, the polymer bulk structure at 700 K and the exclusion grid are used in a Grid-SMD simulation at 700 K to create a cavity (Fig. 4d).
The DW-CNTs are inserted multiple times in the open spaces to render a concentration of ~10% w/w; which corresponds to two DW-CNTs in the PS melt and six DW-CNTs in the PEEK melt.Extra care must be taken if the exclusion grids are in close proximity, as the Grid-SMD simulations may end up opening up bridges of empty space connecting the cavities.In particular, PEEK was a difficult system, not only during the periodic melt creation which required a contraction step, but also during the CNT insertion.For PEEK+CNT system, the CNT insertion has to be gradual: first one DW-CNT at the center of the melt, then three DW-CNTs surrounding the initial DW-CNT, and finally two DW-CNTs at the borders of the periodic cell.After CNT insertion, the PNC structure is heated at 700 K and annealed with cooling steps of −100 K.For annealing of PNCs, the cooling steps were changed from the −100 K/2 ns rate used for melts to a slow-quenching until density convergence at that temperature.This allows a better packing while the polymer is fluid at high temperatures, resulting in a shorter equilibration at 300 K (Fig. 3c, d), and smaller simulation time over the entire annealing cycle (Supplementary Table 1).The Fig. 2 Building polymeric bulk.a Workflow to built a periodic polymeric melt.MD simulation steps are shown in gray boxes.b A PEEK polymer structure is replicated into a 3-D grid with random orientations, each chain is presented with a different graytone.Then, polymer chains are moved toward the grid center using a geometrical collapse procedure, which preserves a distance (1 nm) of empty space among polymers (see Supplementary Methods 2).c Polymer aggregate after geometrical collapse step.Scale bar of 100 nm applies to (b) and (c).Note the volume reduction from the grid order in (b) to the geometrical-collapse aggregate in (c).After that, the polymer chains are force-driven toward the center using a steering MD simulation (force-collapse box in workflow scheme), the resulting polymer melt is placed in a periodic box.d Structure after force-collapse step.The following steps (e-g) are carried out using periodic boundary conditions, the periodic boxes shown as black squares.An MD equilibration at high temperature (equilibration box in workflow) does not result in a homogeneous density (Supplementary Table 1).The long polymers are highly interwined, constraining polymer mobility to fill the empty spaces within the periodic cell.Thus, an expand/contract step is required to reach homogeneous density.e MD expansion, note the corners of the periodic box now filled with polymer.f MD contraction, note homogeneous density and polymers crossing the periodic box.Polymers chains are colored in same graytone in (e) and (f).In the end, an annealing MD cycle is used to reach room temperature.g Final periodic bulk structure.Scale bar of 10 nm applies to (d-g).
final densities for PS/CNT and PEEK/CNT nanocomposites were 1.047 and 1.241 g cm −3 , respectively.

Crystal formation
During the heating stage at 700 K, the polymers are fluid and fill all space around and inside the CNTs.As the temperature decreases throughout the annealing cycles, the polymers transitioned from fluid to solid, resulting in structured phases inside and outside the DW-CNTs.The molecular ordering at the outer convex surface of CNTs has been widely reported for polymers with aromatic rings; such as polystyrene 8 , polyurethane 10 , aromatic polyamides 11 , and DNA 40 .In this case, the polymers are straightened to maximize the contacts between the polymer benzene rings and the CNT surface.In addition, we also observed polymer crystallization inside the DW-CNTs (Fig. 5 and Supplementary Movie 1).For PS, the benzene rings act as pendant groups branching from the main chain, and they can either stack among themselves or with the concave CNT surface.In one case, we obtained crystalline PS segment that closely resembles an α-syndiotactic PS crystal 41 (top half in Fig. 5a): the benzene rings are orthogonal to the CNT axis and oriented towards the center, whereas the main chains are exposed to the CNT wall.Such ordering is comparable to the double helix structure of DNA that is maintained inside CNTs 42,43 .Base stacking between PS and the CNT surface distorts any PS linearity, creating semi-crystalline (top half in Fig. 5b) and amorphous (bottom   1. halves in Fig. 5a, b) regions.For PEEK, there are no pendant groups; moreover, each PEEK monomer contains three aromatic rings and a carboxyl group that stiffen the main chain.Inside the CNT, the PEEK benzenes are in contact with the wall, creating helical patterns (Fig. 5c-h).Similar helicity inside CNTs have been reported by Kumar and coworkers for linear bead-polymers 44 and by Lv and coworkers for all-atom models of polyacetylene, polyparaphenylene vinylene, and polypyrrole 9 .These polymers have a topology similar to PEEK; that is, a predominance of sp 2 bonds that impose planarity in the main chain and the absence of bulky pendant groups.It is worth mentioning that in Lv's study 9 , crystal ordering for PS was not reported, as the PS chains were isolated and not part of a dense polymer melt.
Crystallization inside CNTs has previously been observed using High-Resolution Transmission Electron Microscopy for molecules that cluster together through strong intermolecular interactions; such as inorganic crystals [45][46][47] and small organic molecules with highly conjugated moieties 48,49 .Although the encapsulation of PS inside CNTs has also been visualized by High-Resolution Transmission Electron Microscopy 50 , the crystalline PS arrangement has not been reported yet.We hypothesize that this is due to a current limitation of the experimental method, as electron beams transfer kinetic energy to the flexible polymers 51 , displacing any polymer ordering from equilibrium positions into a melted chain.Overall, these results show that CNT confinement induces ordering for carbon-based polymers, the chain alignments being dependent on the particular topology of each encapsulated polymer; namely the planarity of the backbone and the pendant groups.
Mechanical response of PS+CNT systems Steered Molecular Dynamics (SMD) simulations are used to characterize the mechanical properties of the PS melt and the PS+CNT system.Two layers of materials are selected, one layer is restrained, the other one is pulled away, and the polymers and DW-CNTs located in between are rearranged as the simulation progresses.To deduce relative tendencies independent of the actual velocities, the pulling velocities were varied from 0.1 to 100 nm ns −1 .Fig. 6 summarizes the results.We must underscore that for the PS melt, the pulling layer contains only PS atoms; while in the PNC systems, the pulling layer comprises PS and CNT atoms (see Methods).Such setup was chosen because all systems are torn near the pulling layer, similar to the fracture mechanism of rectangular inorganic tablets without indentations 52 .If CNTs are not included in the pulling layer, the CNTs will remain anchored at the large polymer block during fracture, occluding the possibility to observe reinforcement effects and bridging.
Figure 7 shows the stress-strain (σ − ε) curve for a PS melt using a pulling velocity of 0.1 nm ns −1 .The σ-ε curve reveals a ductile material, characterized by a linear deformation up to a yield point around 0.02-0.03ε, followed by a non-linear strain-hardening region that features inelastic deformations.From the σ-ε curves, we computed Young's modulus (E 0.02 ), stress at fracture (ε f ), strain at fracture (σ f ), and fracture energy (U f ).Young's modulus was calculated within 0.02 ε; that is, in the region with a linear elastic deformation.Figs.6a-d show the results for PS melt (blue) and PS +CNT systems with one (green) and two (red) DW-CNTs in the pulling layer.The systems were pulled along the X-axis (squares), Y-axis (triangles), and Z-axis (circles).Figure 6a shows that adding CNT increases E 0.02 independent of the pulling velocity and pulling direction.Conversely, the properties at fracture, ε f and σ f , do not markedly increase when going from the PS melt (blue) to the PS+CNT systems (green and red).In the case of ε f (Fig. 6b), the data are scattered but reveals the same material extension at fracture point, regardless of CNT reinforcements.In the case of σ f (Fig. 6c), PS+CNT systems shows a slight increment compared to the PS melt.Lastly, the fracture energy U f (Fig. 6d) shows a small increase for the PNC with two DW-CNTs in the pulling layer (red), but overall the values are also widely scatter for different velocities and systems.
Current PNC materials are manufactured using a variety of protocols which result in different CNT features; such as radii, multiwall arrangements and dispersion distribution.Thus, the reported mechanical properties are widely dispersed, from enhancing to even decreasing.The most frequently reported mechanical parameter is Young's modulus; Qian and coworkers 28 have found that PS+CNT of 1% w/w CNT enhances Young's modulus about 36-45%, which is close to our enhancement but with a different CNT concentration.On the other hand, it has been reported that composites exhibit a decrease of Young's modulus at low nanotube loadings 53 .For comparison, PS and PS+CNT specimen were manufactured by two methods: Solvent-Casting and Melt-Mixing.Fig. 8a, b and e, f shows snapshots of the PS and PS+CNT composites.The samples were mechanically tested, confirming that the mechanical properties strongly depend on the sample preparation (Table 1).For example: for Solvent-Casting samples, Young's modulus (E) decreases upon the addition of CNT; whereas for Melt-Mixing samples, it slightly increases.Due to this variability, a direct comparison of the overall magnitudes from experiments and simulations is not straightforward.Instead, we explore the variation of the mechanical response with different loading velocities.We considered the simulation results for PS+CNT systems with 2 CNTs in the pulling layer at 1, 10, and 100 m s −1 (Supplementary Table 2-MD simulations SCc2z1, SCc2z2, SCc2z3) together with the mechanical tests for PS+CNT composites manufactured by SC method at 1.33 × 10 −3 and 3.3 × 10 −4 m s −1 (Table 1-2nd and 4th rows).It is worth noting that the difference in the order of magnitude of experiments and simulations; therefore, we can only draw a qualitative comparison.
As a general scaling law, it is assumed a power-law approximation based on Dynamic Quantized Fracture Mechanics 54 : where P is the mechanical property, v is the pulling velocity, α is an exponent for the scaling law, which is positive if the considered property increases with the velocity and negative if it decreases, and k is a constant with fractional physical units.For Young's modulus (E), α is 0.27 with a square of the coefficient of correlation (R 2 ) of 0.96 (Fig. 9a), suggesting that a scaling law can result in good predictions for E across pulling velocities.Similarly for σ f , α = 0.37 with R 2 = 0.97 (Fig. 9b).Dynamic Quantized Fracture Mechanics predicts a limiting exponent α of 0.5 55 ; thus, E and σ f are quantities that are in line with Dynamic Quantized Fracture Mechanics theory.Conversely, for ε f (Fig. 9c) and U f (Fig. 9d), R 2 values are 0.04 and 0.35, respectively (α = −0.02for ε f and α = 0.23 for U f ).ε f and U f have large variations not only across experimental tests (see errors in Table 1) but also in MD simulated pullings (Fig. 6b, d); showing that for those two mechanical properties, their variations with pulling velocities are rather unpredictable in both experiments and numerical simulations.That is, the fluctuations in both experiments and simulations reflect the stochastic nature of ε f and U f during the rupture process.
Finally, we describe atomic events of the rupture process using a pulling velocity of 1 nm ns −1 .Figure 6e shows an SMD simulation of PS.The PS melt shows a minor extension prior to fracture (ε f ) compared to the initial system (ε 0 ); then, it undergoes rupture by a sudden expansion with only a few polymer chains Fig. 7 Stress-strain curve.Plot shows a typical stress-strain curve.The system is PS under 0.1 nm ns −1 pulling velocity (simulation Sz4-Supplementary Table 2).Abscissa axis is strain (ε).Left ordinate is stress (σ), the values are shown in the oscillating line.X symbol marks the fracture point.Gray area is the fracture energy.Right ordinate axis shows Young's modulus (E), the average values over segments of 0.02 strain are shown as dots joined by straight lines.2).
bridging the fracture zone.Figure 6f, g shows SMD simulations for PNCs with one and two DW-CNTs in the pulling layer, respectively.Both PNC systems show that the presence of embedded DW-CNTs increases the number of polymer chains bridging the rupture zone after fracture.This network of bridging chains is particularly crowded around the DW-CNTs (right panels of Fig. 6e-g).Prior to fracture, the DW-CNTs cannot change orientations within the polymer matrix; whereas after fracture, the open space allows them to rearrange and align along the pulling direction (right column Fig. 6f).Moreover, our SEM images of experimental samples show isolated CNTs and CNT bundles oriented orthogonal to the fracture surface (Fig. 8d, h).Similar alignment has also been reported in experiments 28 .We include a movie of PEEK+CNT pulling (Supplementary Movie 2), where a DW-CNT also aligns perpendicular to the fracture plane.
In this work, we propose a MD simulation workflow based on annealing cycles and grid-SMD simulations to assembly nanocomposites.The protocol involves collapsing a polymer melt into a periodic box, opening volumes for insertion of CNTs, and a final MD annealing to reach constant volume.Our workflow successfully addressed the challenge to build a representative volume element of nanocomposite, which includes CNT reinforcements embedded in PEEK and PS chains of high molecular weight.To ease its usage, we have included the basic operations in a TCLlibrary that can be used with the software VMD.
To summarize, the computational assembly of PS+CNT and PEEK+CNT nanocomposites revealed polymer ordering and crystallization within the CNTs.The polymer topology; namely the planarity of the backbone and the pendant groups, determines the ordering inside the CNTs.Moreover, the PS+CNT models were used to measure Young's modulus (E), strain at fracture (ε f ), stress at fracture (σ f ), and fracture energy (U f ); showing an increase in E due to the addition of CNT reinforcements.Experimental tests on PS+CNT nanocomposites confirm that the mechanical properties strongly depending on the manufacturing process.Thus, the enhancement we observe for our simulation systems suggests that highly dispersed CNTs can in principle strengthen a polymer material, whereas under experimental conditions such ideal dispersion is difficult to achieve.A scaling law derived from Dynamic Quantized Fracture Mechanics can describe the combined experimental and numerical data well for E and σ f .Given its broad applicability; it is expected that Dynamic Quantized Fracture Mechanics can be used to predict the mechanical properties of other polymer nanocomposites of interest.
As expected, the MD workflow is not free of caveats and there is plenty of room for improvement.The high temperatures used MD annealing can modify the stereochemistry of carbon atoms; thus, if the polymer has chiral centers, it must be revised for stereochemical correctness 33 .For example, in Supplementary Fig. 1, we include snapshots of a biological composite: a spider silk model that combines crystalline beta dices and amorphous regions of protein, which was stereochemically corrected during annealing.Another possible application is to use the opening step to build models of spongy polymers and aerogels 56 , where cavities are distributed throughout the material.We also recommend to use only up to 700 K for heating cycles, which provides enough mobility for a carbon-based polymer and avoids atomic collisions that require multiple MD restarts.The MD workflow is still computationally expensive and requires several intermediate steps to reach a well-equilibrated system.Some of  these limitations can be addressed by combining all-atom and coarse-grained simulations 20 , which can also include explicit π-π interactions in the force field.We believe that these and future studies can help to move towards a quantitative prediction and rational tailoring of nano-composite mechanics.

Workflow routines
All the workflow steps are included in a computational tool called "nanocomposite", which is written in TCL language and can be loaded in VMD with the command: package require nanocomposite; the different routines can be called as follows: nnc command [-options] Currently, the "nanocomposite" tool provides the following commands: randomGrid, to create a grid of randomly oriented structures; geoCollapse, to apply the geometric-collapse routine to the grid; and phantomVolume, to map an exclusion grid.Since equilibrium and non-equilibrium MD simulations are used at several building stages (gray boxes in Figs. 2 and 4), "nanocomposite" also includes namdConfiguration, this command eases the creation of configuration files.In Supplementary Methods 1 we provide further details about the "nanocomposite" installation and options as well as a reference to an online tutorial.

Polymeric melt models
The bonding topologies and force field parameters for PS, PEEK, and CNT were obtained by homology from the CGenFF force field 57 .CGenFF is a generic force field that reproduces properties of carbon-based polymers.CGenFF was initially parametrized for organic model compounds to be used together with CHARMM biomolecular simulations.Its key feature is transferability; that is, the model compounds are used as pieces to build more complex materials; thus, it has also been used off-the-shelf to model a variety of polymers; such as poly(ethylene glycol) diacrylate 58 , poly(lacticco-glycolic acid) 59 , poly(N-acryloylglycinamide) 60 , and polyethylene terephthalate 21 .CGenFF has also been shown to reproduce Young's Modulus of graphene materials 61 .NAMD is used for MD simulations, as it provides a parallel implementation for two user-defined force features, namely TclBC 62 and Grid-SMD 38 , which are used to pack polymer melts and insert CNTs, respectively.Details about the MD simulation settings, topology and parameter files are available in Supplementary Notes 1 and 2.
Ten PEEK polymers and ten PS polymers were build using PSFGEN, each polymer chain composed of 100 monomers.Dangling atoms at the polymer ends were methylated.To disorder the linear structures, each polymer was minimized and equilibrated for 10 ns at 1000 K using non-periodic conditions.The polymers were randomly chosen and arranged into a grid of 5 × 5 × 5.After the geometric-collapse step, the blob arrangement was collapsed by a force-driven MD simulation as follows: first, a box is defined at the center of the polymer blob with lengths of 7 nm for PEEK and 5 nm for PS; then, a force of 5 pN towards the blob center is applied to all carbon atoms outside the box.The forcecollapse step was carried out at 1000 K and non-periodic condition (Supplementary Table 1).
The aggregated polymers are inserted into a periodic box of 17.0 nm side for PS and 22.4 nm for PEEK (Fig. 2d) and MD equilibrated at 1000 K in NPT conditions up to filling all empty spaces with polymer.To redistribute the mass at the boundaries of the periodic box, we carried out an expansion/ contraction step (Figs.2e, f); which performs MD simulations with either an outward or an inward force towards the cell center.For the expansion step; an inner cube of side L expand is defined at the center of the periodic cell with L PBC is the length of the periodic box and L free is 1 nm.An outward force (f out ) is applied to all carbons inside the inner cube: where f e is 2 pN and d c2c is the vector distance from the carbon atom to the box center.f out acts stronger at the borders where there is empty space and softer at the center where atoms are crowded.The parameter L free specifies a region near the cell borders where f out is zero; to prevent that f out flips direction when atoms cross the periodic cell.The expansion simulations were carried out at constant volume (NVT) to prevent enlargement of the cell size due to the added forces.For the contraction step, a spherical region of radius 0.35L PBC is defined at the center; then, an radially inward force of 1 pN is applied to the carbon atoms located inside that region.The contraction simulation is carried out at constant pressure (NPT) to allow shrinking of the periodic cell.Once uniform mass Fig. 9 Dynamic quantized fracture mechanics modeling of mechanical properties.Plots show fitting procedures for combined experimental (blue triangles) and simulated (red circles) data.Data points are linearly fitted to a Dynamic Quantized Fracture Mechanics Model (Eq.( 1)) using a log-log scale.E is Young's modulus (a).σ f (b), ε f (c), and U f (d) refer to strain, stress, and fracture energy at breaking point, respectively.
distribution is achieved, annealing cycles are applied up to reaching density convergence at 300 K.

PNC models
In our models, we considered the principal features that characterize a generic experimental sample (multi-walled CNTs, amorphous polymers), while maintaining model sizes of several hundred-thousand atoms and a total MD simulation time of about 1 μs.Nevertheless, we still had to make some simplifications regarding the system size.Thus, we employed the smallest multilayer CNT, which has been reported by Holt and coworkers 63 : a doublewalled CNT with 1 nm inner radius.Furthermore, our efforts are focused on the CNT-polymer interactions; thus, in order to maximize CNT-polymer contact surface and interactions, the DW-CNTs were distributed separately from each other and CNT agglomeration was not considered.
A double-walled CNT (DW-CNT) was created with the VMD-plugin NANOTUBE, with an axial length of 10 nm and chiral indices of (15,15)  and (21,21) with radii of 1.03 and 1.41 nm, respectively.Dangling carbons at the CNT openings were capped with hydrogens.The DW-CNT was minimized and equilibrated for 10 ns at 300 K using non-periodic conditions.The equilibrated DW-CNT structure was inserted multiple times with different orientations into the polymeric bulks to render a concentration of ~10% w/w.
The volume around 1.1 nm of the DW-CNT structure is selected to delineate a exclusion grid (Fig. 4b).A cartesian 3D grid with 0.2 nm spacing is assembled within that volume.In MDFF, the value of the grid nodes (V grid ) are proportional to the mass density and are used as an attractive potential to steer the MD models into the cryo-EM region.To construct an exclusion grid that repels all atoms, we assigned radial values for the nodes; that is, V grid decreases at a rate of 2 kcal mol −1 per nm of distance from the center.The exclusion grid applies a force (f grid ) on all carbon atoms as follows: where w is a scaling factor of 0.05 for PEEK and 0.10 for PS.
The DW-CNT and the polymer melt with cavities are joined into a single system (Fig. 4e) and MD simulated at 700 K until all empty space is filled with polymer, the DW-CNT being positionally restrained with an harmonic constant of 1 kcal mol −1 Å −2 .After all DW-CNTs are inserted, the restraints are removed from the DW-CNTs and the PNC structure is heated at 700 K and annealed up to reaching density convergence at 300 K (Supplementary Table 1).

SMD simulations
SMD simulations require an empty volume along one cartesian axis; thus, the system has free space to stretch without colliding with its periodic image.Specifically, the periodicity in the pulling axis has to be first removed and then extended.For long polymer chains, the removal of the periodicity also shifts the polymers that cross the periodic cell; resulting in regions of low density at the borders.To avoid this problem; we replicate all polymers three times, one copy remains at its original position and the other two are displaced ± one unit cell distance along the pulling axis; then, the polymers located outside the region between the two layers are deleted and the periodicity is restored with an enlarged value.
This procedure of adding empty space along two axes and applying steering forces to the outer layers has been used to study the mechanical properties for a variety of systems, such as silica 64 , aragonite 52 , and graphene 65 .Canonical ensemble (NVT) is used to avoid a barostat control that can reduce the empty space added.The advantage of using SMD to calculate stress-strain curves is that it allows us to observe the fracture mechanisms at molecular level, revealing the initial crack and subsequent CNT re-aligment after fracture.
The restrained and pulling layers are rectangular volumes of 1 nm thick located at the boundaries of the original unit cell.For PS melts, three systems were build, each system allows for pulling along one cartesian axis (X-, Y-, Z-).For PS+CNT systems, the location of pulling layer is shifted inwards to include DW-CNTs in the pulling region.Four systems were build: three of them contain one DW-CNT (3 nm, 4.5 nm, and 3 nm shifts in X-, Y-, and Z-, respectively); the fourth PS+CNT system contains two DW-CNTs (4 nm shift in Z-).The restrained layer is anchored to its initial position through a harmonic potential with a constant of 1.195 kcal mol −1 Å −2 , whereas the pulling layer is steered at constant velocity using a harmonic constant of 1.195 kcal mol −1 Å −2 .As the polymer melts are amorphous and the CNTs are randomly oriented, each cartesian direction exposes different angular orientations for CNTs and different amorphous polymers.
Moreover, the pulling layers contain different CNTs.Thus, each pulling simulation is using a different molecular arrangement, allowing averaging over three simulations when deducing mechanical properties.
The magnitude of the applied force (f SMD ) and the distance between layer centers (L r2p ) are recorded to obtain stress-strain (σ − ε) curves as follows: where A is the cross-sectional area and L 0 is L r2p at time zero.The material fracture is characterized by a sudden drop in the curve f SMD vs time.From the σ-ε curve (Fig. 7), we computed the fracture energy as the area under the curve at the fracture point (U f ), the stress at fracture (σ f ), and the strain at fracture (ε f ).Young's modulus (E) were calculated as a linear slopes for 0.02 ε segments.All SMD simulations were carried out in NVT ensemble and 300 K. Supplementary Table 2 summarizes the SMD simulations.

Experimental samples and mechanical tests
A solution-based process was adopted to prepare the PS (Styron TM 485) composites.Styron TM 485 is an easy flowing PS, offering good impact strength and flexibility, its density is 1.05 g cm −3 .As-received multi-walled CNTs (Nanocyl S.A., trade name Nanocyl NC7000, average diameter 9.5 nm, average length 1.5 μm) were dispersed in chloroform and sonicated for 2 h in a water bath at 20 ∘ C. The nanofiller dispersion was added to the PS solution (0.1 g cm −3 ) in the same solvent in a ratio yielding 2% w/w CNT content (larger CNT concetrations result in aggregation and scarce mechanical properties).Then, the mixture was stirred for 3 min.at room temperature.The obtained solutions were cast on circular Al molds.The composite material was obtained by evaporation of the solvent at room temperature.The mechanical properties of neat PS and PS+CNT composites were measured by a universal tensile testing machine (Lloyd Instr.LR30K) with a 500 N cell at room temperature.The samples were cut into strips of 10 mm × 60 mm × 0.8 mm (Fig. 8a, b).The extension rates were 5 mm min −1 and 20 mm min −1 , whereas the gauge length was 23.5 ± 0.6 mm.Young's modulus has been calculated from the slope in the stress-strain curve between 0.05% and 0.25% strain.Five samples for each composition were tested.Melt mixing and compression molding processes were also used to prepare both PS and PS+CNT composites.As-received PS (Styron TM 637) and CNTs were melt mixed using a Haake Polylab system equipped with an internal mixer at 190 ∘ C and 5.2 rad s −1 for 10 min.Styron TM 637 is a PS with high strength designed for use in extrusion blending, its density is 1.05 g cm −3 .The CNTs were added to yield a 2% w/w proportion.The mixtures were pressed into squared plates of 10 cm side and 2 mm thick using a hot-plate hydraulic press equipped with water cooling system.Then, dogbone shaped samples were cut from the plates (Fig. 8e, f).The mechanical properties of PS and PS+CNT composites were measured by an universal tensile testing machine (Lloyd Instr.LR30K) with a 30 kN load cell at room temperature.The samples were characterized according to the ISO 527-2/5 test method.Five samples for PS and PS+CNT composites were tested.
Field-emission scanning electron microscopy (Zeiss Supra 25) imaging was used to investigate PS+CNT nanocomposite samples.Cross sections were obtained by fracture in liquid nitrogen.
Once loaded, the different routines can be called using the prefix nnc as follows: nnc command [ -options ] Currently, "nanocomposite" provides four commands: • randomGrid : creates a spaced grid of randomly oriented structures.
• geoCollapse : aggregates the grid structures by using translations and rotations.
• phantomVolume : maps an exclusion volume with the shape of a molecule.
Each command provides different options.For the most common tasks, default options are included; for example, during the geometric collapse step (nnc geoCollapse), the aggregated structures end up separated by a default distance of 1 nm.Such distance can be changed using the flag -gapLength to aggregate the grid structures either closer or further apart.Some algorithms and routines within "nanocomposite" can find further applications; for example, the routine emptyLength used for geometric collapse, returns the empty space around a molecule and can be used to quantify the confinement/crowding of a molecule.
To ease the usage of "nanocomposite", a tutorial with the basic examples is available at: https://github.com/nanocomposite/nncTutorial The tutorial material include a PDF guide, and the compressed files nnBuild.zipand nncTest.zipwith exercises to create random grids, perform geometric collapse protocols, generate grids, and set NAMD configuration files.
The "nanocomposite" tool is not limited for polymer-CNT composites and can be used in a myriad of system.For example, following an extended protocol, "nanocomposite" was used to generate three dimensional all-atom structures of two-phase silk protein, which combines amorphous and crystalline phases of the same proteinaceous polymers.
The protocol for creating spider silk starts by creating an amorphous protein block.That is; creating a random grid, reordering the protein polymers with the geometric collapse procedure, and annealing at high temperature.Then, a exclusion volume with the shape of a protein crystal (Supplementary Figure 1.a) is used to create a cavity in the amorphous protein bulk.After that, polymer chains are steered into a crystalline arrangement inside the cavity and the entire structure annealed up to reaching room temperature (Supplementary Figure 1.b).The structure is tested for stereochemical correctness after each annealing step.

Supplementary Notes 1
All MD simulations are performed with NAMD. 2 The MD configuration options are set as follows: 1 fs timestep; short-range interactions (van der Waals and electrostatics) are calculated using a cutoff of 1.2 nm with a switching function starting at 1.0 nm.For periodic systems, long range electrostatics are computed using Particle Mesh Ewald method with a grid density of 1 Å−3 .Minimizations are carried out with a conjugate gradient method implemented in NAMD.MD simulations in NVT ensemble are kept at constant temperature using a Langevin thermostat, whereas simulations in NPT ensemble are kept at 1 atm pressure and a constant temperature using a hybrid Nosé-Hoover Langevin piston.All MD simulations carried out at a temperature higher than 300 K are preceded by a minimization.Supplementary Table 1: Summary of MD simulation campaign to build bulk polymers and polymeric CNT composites.S1-S9, E1-E9, SC01-SC07, and EC01-EC12 labels refer to MD simulations for PS, PEEK, PS+CNT, and PEEK+CNT, respectively.Simulation types (4 th column) are defined in Figures 2 and 4

Supplementary Discussion
Here, we present the structural analysis for the atomic models described in the manuscript (PS, PS+CNT, PEEK, and PEEK+CNT).All systems contains 125 polymer chains composed of 100 monomers each.

PS Bulk
Analysis was performed on a PS bulk model at 300 K (System S9 -Supplementary Table 1).
• Density: We calculated the density of the polymer and also of a half-volume box centered at the origin for the last 1 ns of the trajectory (Supplementary Figure 4).The density of the PS bulk is 1.0023 g cm −3 (standard deviation=0.0007).The density of the half-volume box is 1.0018 g cm −3 (standard deviation=0.0016).Thus, the densities are almost equal for both regions (Confidence level=95%, Student's t-test, H 0 :µ 1 − µ 2 =0, p-value=0.40),revealing that the PS bulk reaches homogeneous density.• Dihedral angle distributions: Dihedral angles were calculated using two atom selections: one takes into account only backbone atoms (Supplementary Figure 5c), the other takes into account backbone and side chain atoms (Supplementary Figure 6c).For both atom selections, the dihedral distributions are similar (Supplementary Figures 5a and 6a).The distributions in Supplementary Figure 5a agree with previous experiments and simulations of amorphous atactic PS. 3 • Radius of gyration: The mean radius of gyration is around 24 Å (Histogram in Supplementary Figure 7).If the polymer chains were streched, half of its length would be approximately 100 Å.Thus, radius of gyration calculation shows that the polymer chains are coiled.The variation of the mean radius of gyration was lower than 0.5 Å throughout the simulation.

Radius of gyration (A °)
Frequency • RMSD: For the RMSD calculation (Supplementary Figure 8) we selected the carbons of the polymer backbone.Supplementary Figure 8: RMSD of the polymer backbone carbons for th system System S9. Supplementary Figure 15: RMSD of the PEEK backbone for the System E9.
• Dihedral angle distribution inside CNTs: For the entire system, the dihedral distribution is similar to the PEEK bulk distribution with a slight decrease in the broad peak centered at 0 • .For the PEEK inside of CNTs, we did not find any significant change in the distributions (Supplementary Figures 12d and 13d).
• Radius of gyration: The mean radius of gyration for all polymer chains is around 54 Å (Histogram in Supplementary Figure 14).
• Radius of gyration inside CNTs: We calculated the radius of gyration inside and outside CNTs.For this calculation, we used nine PEEK chains located inside CNTs, that have at least ten residues inside and ten residues outside the CNTs.The results are listed in Supplementary Table 4.The radius of gyration does not show any specific trend; thus, the PEEK stretching can be similar inside or outside the CNTs.
Supplementary • RMSD: For the RMSD calculation we selected the carbonyl carbons, the ether oxygens and, in the aromatic rings, the carbons that are bonded to oxygens or the carbonyl (Supplementary Figure 17).

Fig. 1
Fig. 1 All-atom polymeric nanocomposite models.The snapshots show PNCs composed of (a) PS+CNT and (b) PEEK+CNT.PS and PEEK polymers are shown in licorice representation, each color represents a different chain.CNTs are shown as cyan beads.The systems are periodic, the borders of the periodic cells being shown as black cubes.To observe the location of CNTs, only half of the polymeric chains are shown.CNT concentration is 10% w/w for both systems, namely two CNTs for PS melt and six CNTs for PEEK melt.A CNT is hidden by PEEK polymers in the bottom back corner.Chemical structures for the polymers are shown next to each MD model.Scale bar of 1 nm applies to both MD systems.

Fig. 4
Fig. 4 Building polymeric nanocomposite.a Workflow to built a PNC.MD steps are show in gray background.Snapshots (b-f ) show insertion of a CNT (b) in a PEEK bulk model.c A 3-D grid with the shape of the CNT structure, so-called exclusion grid.The values in the grid points are radial, the largest value at the grid center.Grid values are shown as spherical isosurfaces.d The exclusion grid is used in a Grid-SMD simulation (opening box in workflow) to create an empty space in the PEEK bulk.e CNT and PEEK bulk structures are combined, followed by equilibration and annealing MD steps.CNT is shown as drak gray beads.f Final PNC structure, note PEEK chains inside the CNT.

Fig. 5
Fig. 5 Polymer ordering inside carbon nanotubes.Snapshots show the polymer arrangement within CNTs for PS (a, b) and PEEK (c-h) melts.Polymers are shown in licorice representation, each color represents a different chain.CNT are presented as cyan beads.Hydrogen atoms are not shown.Scale bar at the bottom right is 1 nm.

Fig. 6
Fig. 6 Mechanical properties of nanocomposites.Plots (a-d) show the mechanical properties for PNC containing PS for different pulling velocities (abscissa axes).PS bulk melt without CNTs is shown in blue.PS+CNT systems with one CNT in the pulling layer are shown in green.PS+CNT system with two CNT in the pulling layer is shown in red.E 0.02 (a) is Young's modulus at 0.02 strain.ε f (b), σ f (c), and U f (d) are strain, stress and fracture energy at breaking point, respectively.Squares, triangles and circles refer to pullings along the X-, Y-, and Z-axes, respectively.Snapshots e-g show MD pulling simulations for PS (e), PS+CNT grabbing one CNT (f), and PS+CNT grabbing two CNTs (g).Polymer chains are presented as gray lines.Restrained and pulling layers are colored in dark gray.CNTs are shown as cylinder.Left, middle, and right snapshots show initial conformations (ε = 0), extensions at breaking point (ε f ), and over-extensions (ε ~1), respectively.Snapshots were taken from pulling MD simulations (a) Sz3, (b) SCc1z3, and (c) SCc2z3 (Supplementary Table2).

Fig. 8
Fig. 8 Experimental samples for PS+CNT nanocomposites.Figures show macroscopic samples for PS (a, e) and PS+CNT (b, f) prepared using Solvent-Casting (a, b) and Melt-Mixing (e, f) manufacturing processes.SEM images show the cross section of the PS+CNT samples prepared using Solving-Casting (c, d) and Melt-Mixing (g, h).For Solving-Casting samples, SEM image at high magnification (d) shows individual CNTs (green arrow) as well as wrapped CNTs with PS (blue arrow).For Melt-Mixing, SEM image at high magnification shows the presence of large bundles of individual CNTs (upper half of panel h).

Supplementary Figure 1 :
Spider silk model.The command phantomVolume creates an exclusion grid with the shape of a molecule.(a) All-atom beta sheet dice (arrows) and the corresponding grid (green isosurfaces).Such grid opens an empty space in a bulk material.(b) Coiled proteins chains (red) are drivenback into the open volume, but the structures are steered towards beta sheet arrangement (yellow).A fitting algorithm chooses the chain orientation, either parallel or antiparallel, based on proximity to the cavity.One protein is highlighted in blue, note the protein gets in and out of the beta-sheet.

Supplementary Figure 4 :
Densities for System S9. Top left: Calculated density for the entire periodic box.Top right: Calculated density for a half-volume box centered at the origin.Bottom: Superimposed plots.

Supplementary Figure 7 :
Histograms of the radius of gyrations for the last frame of PS bulk and the last frame of PS nanocomposite.The plots show 25 bins.

Table 1 .
Mechanical properties of PS and PS+CNT composites prepared by solvent casting (SC) and melt mixing (MM) for different pulling velocities.
Young's modulus (E), stress at breaking point (σ f ), strain at breaking point (ε f ), and fracture energy (U f ) are listed in 5th, 6th, 7th, and 8th columns, respectively.a Manufacturing process.b CNT concentration in % w/w.c Pulling velocity in mm min −1 .
in the manuscript.Cell sizes (9 th column) are reported for the last frame of the MD simulation.Abbreviations Geo.Col.(2nd column), Force Coll.(4 th column), Equil.(4thcolumn), gForce (5 th column), and Non.Per.(6thcolumn) refer to Geometric Collapse, Force Collapse, Equilibration, Grid-SMD, and Non Periodic, respectively.Systems did not reach homogeneous density and were discarded (see workflow in Figure2in the manuscript). a

Table 4 :
Radius of gyration for 10 residues of the same polymer chain inside and outside the CNT.Data is of the last frame of the PEEK NNC trajectory.Not all polymer chains that have atoms inside the CNTs appear in this table.