Investigation of the Impact of Cross-Polymerization on the Structural and Frictional Properties of Alkylsilane Monolayers Using Molecular Simulation

Cross-linked chemisorbed n-alkylsilane (CH3(CH2)n−1Si(OH)3) monolayers on amorphous silica surfaces have been studied and their structural properties and frictional performance were compared to those of equivalent monolayers without cross-linkages. The simulations isolated for the first time the effects of both siloxane cross-linkages and the fraction of chains chemisorbed to the surface, providing insight into a longstanding fundamental question in the literature regarding molecular-level structure. The results demonstrate that both cross-linkages and the fraction of chemisorbed chains affect monolayer structure in small but measurable ways, particularly for monolayers constructed from short chains; however, these changes do not appear to have a significant impact on frictional performance.


Molecular Dynamics Simulations
Molecular dynamics simulations of fully and partially chemisorbed alkylsilane monolayers featuring cross-linkages have been performed under equilibrium and nonequilibrium conditions. Simulations were conducted using the optimized potentials for liquid simulations all-atom (OPLS-AA) force field [68]. Parameters used for this work were taken from Lorenz et al. [52] for silica and Jorgensen et al. [68] for alkanes. Some minor adjustments were made to these parameters, as detailed in Tables S1 and S2.  Simulations of chemisorbed alkylsilane monolayers were performed with the LAMMPS simulation engine [70] using a time step of 0.5 fs. Temperature was controlled via the Nosé-Hoover thermostat [71,72] with a temperature damping parameter of 50 ps. The equations of motion were integrated using the multiple time step algorithm rRESPA with time steps of 0.3 fs for bonds, 0.6 fs for angles (valence and dihedral), and 1.2 fs for Lennard-Jones and electrostatic interactions, which were computed using the particle-particle, particle-mesh (PPPM) algorithm for slabs (i.e., electrostatic interactions were not calculated across the nonperiodic boundary). Lennard-Jones interactions were computed using a cutoff radius of 10 Å, in accordance with prior simulation studies [14,16,23,24,58,59]. The positions of atoms in the silica surfaces were not integrated through time, so the relative positions of atoms within a substrate remained constant (i.e., the silica substrates behaved as rigid bodies). For the nonequilibrium simulations, monolayer-coated surfaces were first mirrored over the surface plane (i.e., the xy-plane) to create two opposing monolayers, which were then compressed at a rate of 1 m/s until coming into contact. Four simulation snapshots were taken from each compression trajectory and used as starting configurations for independent shearing runs at fixed separation distances. Constant velocities of +5 and -5 m/s were then applied to the upper and lower silica surfaces, respectively, in the x-direction. The friction force (i.e., sum of all forces in the xdirection) and normal load (i.e., sum of all forces in the z-direction) on each monolayer-coated surface were determined periodically over time.
Partially chemisorbed alkylsilane monolayers were parameterized using the Foyer atom-typing package [73], and simulations were performed using the Gromacs molecular dynamics engine (v5.1.0) [74]. Hydrogen bonds were held fixed using the LINCS algorithm [75], removing high frequency atomic motions, allowing for a time step of 2.0 fs to be utilized. Temperature was controlled via the Nosé-Hoover thermostat [71,72] with a time constant of 1.0 ps. Equations of motion were calculated using a leap-frog integrator [76], and electrostatic interactions were computed using the particle-mesh Ewald technique [77,78] with a correction for slab geometries. For equilibrium simulations, all atoms were allowed to move except for those in the outer 0.4 nm of each surface; however, for shear simulations, atoms in the bottom surface were only allowed to move in the z-direction (normal to the monolayer plane), with a constant normal force applied to these atoms. Constant normal loads of 5, 15, and 25 nN were considered. To induce shear, atoms in the top surface were allowed to move only in the x-direction (the direction of shear) and this surface was coupled via a force constant of 240 kcal/mol/Å2 to a fixed point moving at a velocity of 10 m/s, yielding shear at a constant average velocity. Flow control of simulations and analysis of partially chemisorbed monolayers was achieved using the Signac and Signac Flow packages [79]. The full MoSDeF workflow for these simulations, including initialization, analysis, and run scripts, is hosted on Github and is available at https://github.com/summeraz/xlink-screening.
We note that the simulations of chemisorbed monolayer systems performed at constant separation yield standard deviations in normal load below 5%, so under the conditions of this study, "constant separation" and "constant normal load" ensembles appear to be mostly equivalent.

Determination of System Properties
For monolayer systems at equilibrium and undergoing shear, the nematic order parameter (S2) and average tilt angle (q) were calculated using the following approach. First, average tilt was determined from the characteristic vector, û, describing the long axes of the alkylsilane monomers; ûk for monomer k can be found from the moment of inertia tensor of k, Ik, where each component of the 3 × 3 matrix is given by "#,% = ∑0)12 )* ), #% − )# )%/. (S1) N is the total number of atoms in the monomer, mi is the mass of atom i, ri is the position of i relative to the molecule's center-of-mass, and dab is the Kronecker delta. ûk is the eigenvector associated with the smallest eigenvalue of Ik [80]. Given û, S2 can be calculated by constructing the second-rank ordering tensor for each monomer k, Qk, where each element of the 3 × 3 matrix is given by The director of molecule k, nk, is given by the eigenvector associated with the largest eigenvalue of Qk. The order parameter S2 is , = 〈 cos , − 〉 where cos qk = nk × ûk, and the brackets 'áñ' denote an average over all monomers [81].
The gauche defect fraction was also calculated for monolayers at equilibrium. In alkylsilane monomers, the dihedral angle (ϕ) for twisting a C-C-C-C quadruplet about the central C-C bond is given by (S3) where b1, b2, and b3 are vectors representing the three bonds in the quadruplet; b2 corresponds to the central bond [83]. A dihedral angle is considered to be in the trans state if ϕ is between 90° and 270°, while angles outside this range are considered to be gauche defects [82]. The gauche defect fraction is defined by the fraction of total C-C-C-C dihedral angles that are in the gauche state.