Molecular modeling and simulation of atactic polystyrene/amorphous silica nanocomposites

The local structure, segmental dynamics, topological analysis of entanglement networks and mechanical properties of atactic polystyrene - amorphous silica nanocomposites are studied via molecular simulations using two interconnected levels of representation: (a) A coarse - grained level. Equilibration at all length scales at this level is achieved via connectivity - altering Monte Carlo simulations. (b) An atomistic level. Initial configurations for atomistic Molecular Dynamics (MD) simulations are obtained by reverse mapping well- equilibrated coarse-grained configurations. By analyzing atomistic MD trajectories, the polymer density profile is found to exhibit layering in the vicinity of the nanoparticle surface. The dynamics of polystyrene (in neat and filled melt systems) is characterized in terms of bond orientation. Well-equilibrated coarse-grained long-chain configurations are reduced to entanglement networks via topological analysis with the CReTA algorithm. Atomistic simulation results for the mechanical properties are compared to the experimental measurements and other computational works.


Introduction
The complexity of intermolecular interactions and confinement in polymer-nanoparticle systems leads to spatial variations in structure and dynamics at both the meso-and nanoscale, which are of the order of 100 nm -1 µm and 1 -100 nm, respectively. Nanomaterials fabricated by dispersing nanoparticles in polymer melts have the potential for unique optical, electrical, rheological, and permeability properties, as well as mechanical performance that far exceeds that of traditional composites [1,2,3,4]. The attributes of a heterogeneous material depend on the properties of its component materials and on the characteristic length scale of their coexistence, as well as on the chemical and morphological details of the dispersion. Measurements of the rheological properties of polymer matrix nanocomposites in the melt have revealed unexpected phenomena of great practical importance for the molding technology of polymers. Of particular interest are the experiments of Mackay et al. [1,2,3], who showed that the viscosity of polymer melts with dispersed nanoparticles may be lower than that of their pure counterparts.
Kumar and Krishnamoorti [5] has presented in a review article that despite intense research activity on polymeric matrix nanocomposites, the exact mechanisms by which nanoparticles affect the properties in often unexpected ways still remain largely unknown. The importance of nanocomposites has been known for some time, e.g. in the use of carbon black for modifying the properties of elastomers for car tires. More recently, interest in these materials was stimulated by the discovery [6] that the incorporation of small amounts of nanoparticles in a polymer matrix can result in a simultaneous increase in both the elastic modulus (stiffness), and the energy required for their deformation to failure (toughness), contrary to what happens in the case of filler particles with macroscopic (>1mm) dimensions.
The dynamics of polystyrene nanocomposites has been studied at a microscopic level using neutron scattering techniques. Neutron scattering measurements offer the unique possibility of analyzing the spatial arrangement and dimensions of individual components and provide an excellent tool for assessing the microscopic dynamics of nanocomposites at various length scales. Using quasi -elastic neutron scattering (QENS) measurements, nanocomposites of polystyrene (PS), poly (methyl methacrylate) (PMMA) and tetramethylbisphenol -A polycarbonate (TMPC) with silica and fullerenes have been studied. Nanoparticles were found to bring about an increase in the glass transition temperature of these materials, as measured by differential scanning calorimetry (DSC) and dynamic mechanical spectroscopy (DMA). These authors suggested that the suppression of the local relaxation dynamics of the nanocomposite is consistent with an enhancement of cohesive interactions in the system, which may be the root of the increase in T g for the nanocomposites. A very small increase of T g for silica -polystyrene nanocomposites was observed by Kontou and Anthoulis [7].
Vogiatzis et al. [8] have investigated the structure of a polystyrene matrix filled with tightly crosslinked polystyrene nanoparticles, using Monte Carlo simulations of a coarse -grained model. Their work shows that the local structure of the polymer matrix in the vicinity of the nanoparticle is different in many ways from that of the corresponding bulk, both at the segment and the chain level. The local polymer density profile near the particle displays a maximum and the bonds develop considerable orientation parallel to the nanoparticle surface. Using atomistic MD simulations, Ndoro et al. [9,10] studied the structure and dynamics around silica nanoparticles dispersed in a matrix of low molecular weight polystyrene. They observed that near a nanoparticle a polymer layer of increased density was formed, due to favorable interaction energy. Polymer bonds in this layer were highly oriented parallel to the surface of the particle, and segmental dynamics was slower relative to the bulk. The effect of the nanoparticle curvature and grafting density on the mean square displacement of free polystyrene chains and also on the mean relaxation time of various intramolecular vectors was investigated as a function of separation from the nanoparticle surface. Confinement, reduced surface curvature, and densification resulted in a reduction of the mean square displacement and an increase in the mean relaxation time of the C a -H bond vector and chain end-to-end vector in the vicinity of the surface.
Brown et al. [11] have examined the structure and dynamics of a system containing an inorganic nanoparticle embedded in a polymer matrix. In their work, results are presented regarding the variation of the structure and dynamics with distance from the polymer -nanoparticle interface and as a function of pressure.
Li et al. [12] observed that highly entangled polymer chains were disentangled upon increasing the volume fraction of spherical non -attractive nanoparticles. The authors report a critical volume fraction controlling the crossover from polymer chain entanglements to nanoparticle -induced entanglements. While below this critical volume fraction the polymer chain relaxation accelerates upon filling, above this volume fraction the situation reverses.
Using a quasi-harmonic approximation, Lempesis et al. [13] have calculated and predicted mechanical properties of atactic PS in the elastic regime. The authors report a very good agreement of Young's modulus and Poisson's ratio with experimental measurements [27].

Experimental Methods
The simulation of the nanocomposites is based on a two-tier approach. Equilibrating the melt takes place at the coarse-grained level, where each polymer interaction site corresponds to a monomer unit and each nanoparticle is represented as a single sphere, using powerful connectivity -altering MC algorithms [16]. In particular, a polymer coarse-grained site stands for a PS diad (meso or racemo). PS coarse-grained sites are governed by bonded and non-bonded interactions derived from the atomistic force field via Iterative Boltzmann Inversion (IBI). Silica nanoparticles interact with coarse-grained sites through Hamaker integrated potentials derived directly from the atomistic LJ parameters [17]. Next we restore the atomistic detail through a process which is based on a combination of local MC moves, quasi-Metropolis reconstruction of atomistic sites and gradual minimization of the potential energy of the generated configuration. Specifically, reverse mapping of coarse-grained PS to unitedatom structures is applied retaining implicit coarse-grained tacticity (1 coarse-grained site = 8 united atom carbons) [18].
Atomic-level interactions between polymer chains and nanoparticles play a central role in this work. Given the low polarity of PS segments, they were assumed to consist of excluded volume and London dispersion contributions, represented by Lennard-Jones (LJ) potentials. LJ parameters between polymer atoms and fullerene carbons were obtained using the Lorentz-Berthelot rules. For polymer-silica interactions the parameterization was based on previous work on zeolite/hydrocarbon systems which has yielded excellent predictions for sorption isotherms and diffusivities [19]. Any chemical bonding between nanoparticle surfaces and polymer was disregarded. Silica nanoparticles of prescribed diameter were generated by excising a spherical region from bulk amorphous silica simulated with MD and terminating dangling bonds at the surface as silanol groups, following ref [11].
Finally, to study the dynamics of the nanocomposites, we perform MD simulations in the atomistic representation using a united-atom PS forcefield [20] integrated into the LAMMPS [21] MD application.
The PS nanocomposite systems studied consist of silica particles molecules dispersed in a melt of atactic polystyrene. A single particle of diameter 3 nm, 6 nm or 12 nm, surrounded by polymeric melt composed of chains of molecular weight 20 to 200 kg/mol is simulated under nearly infinite dilution conditions. The coarse-grained MC simulations are performed in the canonical ensemble (NVT) at a temperature of 500 K, for 8 up to 20×10 9 moves. Subsequent reverse mapping to the atomistic level takes 590×10 6 MC moves. The MD simulation starts in the isothermal-isobaric (NPT) statistical ensemble at 500 K and 1 bar and continues in the canonical and microcanonical (NVE). The integration time step in the MD is set to 1 fs (10 -15 s).

Results and Discussion
The estimated root -mean -square radius of gyration 1/2 2 g R as a function of molar mass for neat and nanocomposite polystyrene systems is shown in Fig. 1. Very good agreement is observed for all the molar masses examined for neat polystyrene systems. As far as the nanocomposite polystyrene systems are concerned, the presence of the nanoparticle seems to affect slightly the root -meansquare radius of gyration.
A fusion of trans and gauche conformational states is apparently produced by the atomistic potential at T = 500 Κ. Vogiatzis and Theodorou [18] observed that torsion angle distributions in their reverse -mapped structures are extremely close to those obtained from a 80mer structure that was directly equilibrated by MD at the atomistic level, without intervention of any coarse -graining and reverse mapping. This observation implies that one can achieve well -equilibrated structures, down to the atomistic scale, using the reverse mapping procedure. In the case of racemo diads, the percentage of gauche states increases slightly in the presence of the silica nanoparticle. This effect may be attributed to the effort of the chains to engulf the silica nanoparticle. The overall content of trans conformations and the total percentage of g conformations in the reverse -mapped structures are in reasonable agreement with experimental NMR measurements [23].
Local structure around the silica nanoparticles has been quantified. At short distances from the surface of the nanoparticle, the local density is strongly perturbed. The radial distribution of polymer density around the nanoparticle is seen to exhibit two peaks before stabilizing at its bulk value. Profiles from the atomistic and from the coarse-grained simulation are in good agreement for the same particle diameter, underlining the consistency of the two approaches [24]. Average root mean squared radius of gyration of the coarse-grained chains of neat and nanocomposite polystyrene systems as a function of molar mass in the melt at 500 K (red, green and magenta rhomboid symbols). Neutron scattering measurements for high molar mass polystyrene are also included (blue rhomboid symbols) [22].
Based on the MD trajectories, segmental and local dynamics in the polymer have been characterized. Directly accessible are the orientational decorrelation functions of the phenyl stems. Characteristic times extracted from these decorrelation functions have been compared to dielectric spectroscopy (DS) measurements. In Fig. 2, a slight influence of the presence of the silica nanoparticle (dispersed at infinite dilution) is observed on global segmental dynamics. In this case, the first Legendre polynomial P 1 (t) is employed as a measure of the polymer segmental dynamics. This is a very common measure by calculating time-correlation functions of a vector, v b , representing a chemical bond. This measure of the segmental dynamics of a-PS can be compared with NMR data. It is quite common to fit the long -time behavior of orientational autocorrelation functions of this kind by a modified Kohlrausch -Williams -Watts (mKWW) function [25]. Here, orientational relaxation of a-PS C-CM phenyl ring -backbone vector is measured using P 1 (t) and characteristic times are calculated for bulk PS and PS/SiO 2 systems at volume fraction Φ 1 = 1% and Φ 2 = 6% respectively. We need to extract these characteristic times at different temperatures in order to calculate glass transition temperature Τ g of bulk PS and PS/silica nanocomposites.
The study of local segmental dynamics, when a nanoparticle is present, requires a discretization of space, so as to quantify possible a nanoparticle effects. We choose to carry out an analysis based on a discretization of the simulation box using spherical shells, with the silica nanoparticle acting as the center of the shells. Shells surrounding the nanoparticle are of equal thickness. They were constructed to define different spatial regions as a way to investigate how far the silica surface influences the a-PS segmental dynamics. We clearly observe slower dynamical behavior of this vector near the surface of the silica nanoparticle. The local influence of the presence of the silica nanoparticle on the segmental dynamics of PS are shown in Fig. 3.   Topological analysis of coarse-grained neat and nanocomposite systems of PS has been performed by using the CReTA (contour reduction topological analysis) algorithm [26]. For a neat polystyrene system with 40 chains of 2000 diads each, N e was estimated to be 126 diads, and since the average molar mass of a diad is 104.15 g mol -1 , the corresponding entanglement molar mass M e is 13097 ± 200 g mol -1 . This value is in very good agreement with experimental estimates based on plateau modulus measurements [34]. For a nanocomposite polystyrene system with the same number of chains and molecular weight, the presence of silica nanoparticle with diameter 6 nm increases the N e [24] from 126 diads, which is the value of neat polystyrene system, to 152 diads and the corresponding entanglement molar mass M e accordingly. After calculation of structural and dynamical properties, elastic constants such as the Young's modulus are studied using stress -strain computational experiments at room temperature on bulk PS and PS/SiO 2 systems at volume fraction Φ 1 = 1% and Φ 2 = 6% respectively. For bulk PS, Young's modulus was estimated to be E = 3.7 ± 0.2 GPa. This value is in a good agreement with experimental measurements [27] and other computational works [13]. Likewise, the room temperature Young's modulus for polystyrene -silica nanocomposite systems was calculated from the slope of a line fitted to the stress-strain simulation data at small deformations. The results were E Φ=1% = 4.9 ± 0.3 GPa and E Φ=6% = 5 ± 0.4 GPa, respectively. As we can see, an increase of the Young's modulus relative to bulk PS is brought about by the presence of the silica nanoparticles. Also, the bulk modulus k can be calculated using the strain fluctuations methodology [14,15]. For bulk glassy PS at 300 K, the bulk modulus was estimated as k PS = 3.56 ± 0.4 GPa and for the nanocomposite polystyrene systems as k Φ=1% = 3.75 ± 0.5 GPa and k Φ=6% = 3.94 ± 0.6 GPa. The result for bulk PS is in good agreement with the experimental value k exp,PS = 3 GPa [27].

Conclusion
In this work we have presented results on the local structure, segmental and local dynamics of polystyrene nanocomposites. Developing a hierarchical approach involving two interconnected levels, it was possible to prepare fully equilibrated configurations with high molecular weight chains. The dynamics of the nanocomposite was slower than in the pure melt at a temperature range above the glass transition temperature, in good agreement with available experimental observations. Local density around the nanoparticles was disturbed over a region of thickness of a few nanometers. Entanglements analysis based on the CReTA algorithm shows that the presence of silica nanoparticles increases entanglement molar mass M e , either the silica nanoparticle is treated as a phantom particle or is allowed to interact with reduced chains. The presence of silica nanoparticles in a polystyrene matrix increases Young's modulus value of pure polystyrene. Methodologies for calculating the shear and bulk modulus of our systems through volume and strain fluctuations in the undeformed state have also been used. Our results are in very good agreement with available experimental measurements and other computational works.