Chemically Specific Multiscale Modeling of the Shear-Induced Exfoliation of Clay–Polymer Nanocomposites

We recently showed, using chemically specific modeling and simulation, how the process of intercalation of polymers within clay sheets occurs, transforming the large-scale materials properties by a specific set of spatial and temporal processes that can lead to exfoliation. Here, we use the same hierarchal multiscale modeling scheme to understand the processes that occur during the shear-induced processing of clay–polymer nanocomposites. For both hydrophobic polymers (polyethylene) and hydrophilic polymers (poly(ethylene glycol)), we used free-energy methods to identify the lowest-free-energy separation of the clay sheets; the polymer molecules spontaneously intercalate into the clay interlayer from the surrounding polymer melt. We apply shear forces to investigate exfoliation and find that while exfoliation is promoted by shearing, it is the surfactant molecules that play the dominant role in resisting it.

In this supplementary information, we provide additional discussions of the methods used in our study and present the density profiles within the clay interlayers referred to in the main manuscript.

Simulation protocol
To calculate the free-energy of clay layer separation perpendicular to the basal plane, we created a two clay-layer system (see Figure 1 in the main manuscript) immersed in a box of polymer. The clay sheets are initially placed with the clay layers lying in the xy plane. Each clay layer consists of 672 CG atoms, and has lateral dimensions of approximately 140 × 100 Å 2 . The simulation box was initially 200 × 200 × 200 Å 3 . To examine the role of clay surface charge density, which in turn affects the density of charge-balancing surfactants on the clay surface, we have considered different isomorphic substitution rates in the clay framework.
These correspond to substituting Al 3+ ions with Mg 2+ ions in the octahedral sheet of the clay framework. As described in the main manuscript, we have considered several combinations of surfactants, polymers and substitution percentages; these are listed in Table 1 in the main manuscript. To simulate the free-energy of clay-layer separation, the distance between the z component of the respective centres of mass perpendicular to the clay layers is used as the reaction coordinate during the free-energy calculation. A schematic diagram illustrating the setup is shown in Figure 1.
The procedure to simulate the free energy of clay layer separation is listed below: 1. An un-intercalated clay model was built using an atomic description of the clay layers in vacuum and was subsequently geometry optimised. A short molecular dynamics run (100ps) was then performed, with the clay layers only allowed to move in the z direction. This gives us an initial unintercalated system.
2. The final snapshot from the atomistic simulation was converted into coarse-grained form, and the coarse-grained polymer was subsequently built around the clay platelet using a S2 Monte Carlo procedure to produce a low energy state, drawing on the methods of Theodorou and Suter 1 and the look-ahead procedure developed by Meirovitch 2 . The details are given in our previous publication 3 . This method generates an amorphous polymer system in a relatively low energy state and avoids high energy overlaps. Each polymer was composed of 100 monomer units. Between 599 and 717 polymer molecules were added to the simulation box of 200 × 200 × 200 Å 3 , to give a polymer density of 0.62 g/mL, assuming 20% of the simulation box volume is occupied by the clay platelet. 3. To equilibrate, the clay layers were held fixed during the molecular dynamics simulation and heated up to 500K and with a pressure of 100 atmospheres (to represent melting intercalation conditions) over 60ns.
4. Subsequently, the clay layers were treated as rigid bodies, and only allowed to move in the z direction, thereby ensuring we are only examining the free-energy of clay sheet separation perpendicular to the the basal plane of the clay sheets. 5. A spring was applied between the two rigid layers, initially with a stiff spring constant of 100 kcal mol −1 , increasing by 0.25 Å every 3ns. This effectively fixed the clay-separation, and allowed the polymer and surfactant molecules to equilibrate to the specified clay layer separation. The separation was varied from the un-intercalated value to approximately 60 Å.
6. For each clay-layer separation, the value of the spring was relaxed to 4.7 kcal mol −1 , run for a further 3ns and the resulting separation between the layers recorded over the final nanosecond.
7. The free-energy of clay layer separation is computed using the WHAM procedure 4 . The spring constant used as input to the WHAM procedure was divided by the number of clay CG atoms to give a per CG-atom value for the free energy, (i.e. independent of the number of clay atoms / size of the clay layer. ) All the timescales reported here and in the main manuscript are the explicit timescales from the molecular dynamics simulations; due to the softness of coarse-grained potentials, the dynamics are often much faster. In our previous publication 3 , we estimated this speed-S3 up through a comparison of mean-squared displacement at the atomistic and coarse-grained levels and concluded that for these systems, the speed-up was by a factor of approximately 5.
To simulate the force required to exfoliate one clay sheet from a two-layer tactoid, we used the following procedure: 1. Using the lowest energy separation for each system, we applied a force in the x direction to every atom in the upper-clay sheet. The position of the bottom clay sheet was kept fixed.
The clay sheet was treated a rigid molecule that could only move in the xy plane.

A continuous force is applied at each timestep to the atoms in the upper clay sheet.
Rupture is defined as occurring when the difference in the x component of the centres of mass of the two clay sheets is half the diameter of the clay flakes, i.e. 70 Å; simulations were either continued for 70 ns or until rupture had occurred.

Density profiles
The density profiles perpendicular to the clay surface within each tactoid were calculated at the equilibrium d-spacing determined by the free-energy calculations. The density profiles were calculated using the following procedure: A plane was fitted to the lower clay sheet; the CG coordinates from the simulation were multiplied by the rotation matrix that transforms the fitted plane onto the xy plane. The density profiles were then computed as the number of atoms of a given CG species between planes at z and z + dz, of CG atoms with x and y coordinates that were within the footprint of the clay sheet.