Nanocomposite of Fullerenes and Natural Rubbers: MARTINI Force Field Molecular Dynamics Simulations

The mechanical properties of natural rubber (NR) composites depend on many factors, including the filler loading, filler size, filler dispersion, and filler-rubber interfacial interactions. Thus, NR composites with nano-sized fillers have attracted a great deal of attention for improving properties such as stiffness, chemical resistance, and high wear resistance. Here, a coarse-grained (CG) model based on the MARTINI force field version 2.1 has been developed and deployed for simulations of cis-1,4-polyisoprene (cis-PI). The model shows qualitative and quantitative agreement with the experiments and atomistic simulations. Interestingly, only a 0.5% difference with respect to the experimental result of the glass transition temperature (Tg) of the cis-PI in the melts was observed. In addition, the mechanical and thermodynamical properties of the cis-PI-fullerene(C60) composites were investigated. Coarse-grained molecular dynamics (MD) simulations of cis-PI-C60 composites with varying fullerene concentrations (0–32 parts per hundred of rubber; phr) were performed over 200 microseconds. The structural, mechanical, and thermal properties of the composites were determined. The density, bulk modulus, thermal expansion, heat capacity, and Tg of the NR composites were found to increase with increasing C60 concentration. The presence of C60 resulted in a slight increasing of the end-to-end distance and radius of the gyration of the cis-PI chains. The contribution of C60 and cis-PI interfacial interactions led to an enhancement of the bulk moduli of the composites. This model should be helpful in the investigations and design of effective fillers of NR-C60 composites for improving their properties.


Introduction
Natural rubber (NR), mainly consisting of cis-1,4polyisoprene (cis-PI) with high elasticity, is acknowledged to be virtually irreplaceable in applications such as tires and seals. However, the neat NR usually suffers from poor mechanical strength [1,2]. To alleviate this problem, reinforcing fillers are widely used to improve the mechanical properties the development of CG models for polymers and polymer melts, including the works of Akkermans and Briels, who constructed a CG model with a focus on the preservation of thermodynamic properties [45], followed by a number of studies on the systematic coarse-graining of polymers and melts [46,47], modification of the dissipative particle dynamics (DPD) method for melts [48], and fluctuating soft-sphere models [49] with the aim of improving the accuracy of CG polymer melt modeling. Regarding the interactions of fullerenes and polymers, Huang et al. [50] developed a polymer-fullerene CG model for polymer-based solar cells, and Volgin et al. [51] investigated the diffusion of fullerenes using a CG model. Although there is no lack of CG models, the MARTINI force field [52,53] has become by far the most dominant one, and it has been used for simulations of amino acids, water, phospholipid membranes, fullerenes, some polymers, and RNA [54][55][56][57][58][59][60][61].
Recently, we introduced a MARTINI CG model for cis-PI to study natural rubber [62]. The interactions of the cis-PI and fullerene, however, were not included in the model. That is done in the current work by combining the CG models of cis-PI and fullerene to study cis-PI-C 60 composites. The parameters are based on the MARTINI force field version 2.1 [52]. The thermodynamic and mechanical properties of the cis-PI-C 60 composites were calculated and compered with the experiments and atomistic MD simulations. The influence of the C 60 concentrations on the structural, mechanical, and thermal properties of the composites was analyzed. The results show that the interfacial interactions between the C 60 and cis-PI play an essential role in the composites' mechanical properties.

Molecular Dynamics Simulations
In this study, cis-PI in melts and cis-PI composites at different C 60 concentrations were investigated by coarse-grained (CG) molecular dynamics (CGMD) simulations. The CG model for cis-PI was developed based on the MARTINI force field version 2.1 (University of Groninge, Groningen, The Netherlands) and is described in detail in Table S1. A coarsegrained or superatom bead was used to represent an isoprene monomer. CG mapping of stretching and bending parameters for the cis-PI beads from united-atom MD trajectories were published in our previous work [62]. All systems consisted of 500 chains of 32-mers cis-PI, which is equivalent to 16,000 CG beads in total (an isoprene monomer was mapped into one bead in the CG model).
C 60 molecules were randomly placed into the cis-PI melts with the number of 53, 107, 213, and 427 molecules corresponding to concentrations of 4, 8, 16, and 32 parts per hundred of rubber (phr), respectively. The details of all simulations are shown in Table S2. The CG fullerene model for the C 60 molecule called F16 was taken from Ref. [63]. After steepest descent energy minimization, CGMD simulations were performed using the GROMACS 5.1.1 (University of Groningen, Royal Institute of Technology, Groningen, The Netherlands) package [64]. The number of particles, temperature, and pressure were kept constant (NPT ensemble). cis-PI and C 60 molecules were thermostatted separately at 300 K using the Parrinello-Donadio-Bussi velocity rescale algorithm with a time constant of 30 ps [65,66]. Pressure was held constant at 1 bar by the Parrinello-Rahman barostat [67] with a time constant of 1 ps and compressibility of 4.5 × 10 −5 bar −1 . The reaction-field method was used for the long-range electrostatic interactions [68,69], Lennard-Jones interactions were cut off at 1.1 nm, and periodic boundary conditions were applied in all directions.
The simulation protocol had been tested and used in several prior studies [38,[70][71][72][73]. All production simulations were run with time step of 30 ps for 21 µs. To confirm that equilibration had been reached, time evolution for the end-to-end distance and the radius of gyration were calculated (see Figures S1 and S2). The last 5 µs were used for data analysis and errors were estimated using standard deviation (SD). Molecular visualizations were done using the Visual Molecular Dynamics (VMD) software (University of Illinois Urbana-Champaign, Urbana and Champaign, IL, USA) [74].

Calculation of Thermal and Mechanical Properties
Glass Transition Temperature (T g ) The final states of the simulations after 21 µs were used as the initial structures to calculate the glass transition temperatures (T g ). The systems were cooled from 300 K to 100 K with 10 K interval. The cooling rate was 0.1 K.ns −1 . To equilibrate the systems, additional 500 ns were simulated at each temperature in the NPT ensemble. T g was calculated using the intersection of the two lines describing the behavior at low and high temperatures on density-temperature (ρ-T) curves [75][76][77]. Density versus temperature curves for cis-PI in melts and composites are shown in Figure S3.
The thermal volume expansion coefficient (γ), specific heat capacity (c p ), and bulk modulus (κ) were calculated from the fluctuations of the simulation box volumes (V) and enthalpy (H) using: where k B is the Boltzmann constant, N A Avogadro's number, and m is the total mass of the system in atomic mass units. The angular brackets refer to averaging over simulation time.

Calculation of Solvation Free Energy
Thermodynamic integration (TI) [78] was applied to calculate solvation free energies of C 60 in water, pure cis-PI, and cis-PI-C 60 composites. The last frame from each of the 21 µs simulations was used as the initial structure. A C 60 was added into the system and used as a coupling molecule. After steepest descent energy minimization, the MD simulations with the coupling parameter (λ) between 0 and 1 with interval of 0.05 were performed under the NPT ensemble. All systems were run for 300 ns, and the last 100 ns was used for analysis. The Bennett acceptance ratio (BAR) was applied to estimate the solvation free energy [79]. The solvation free energies are shown in Table S3. The effect of C 60 concentration is discussed in Section 3.4 ( Figure 5).

Coarse-Grained Model Based on MARTINI Force Field of cis-1,4-Polyisoprene
The macroscopic and structural properties of the cis-PI in the melts were analyzed and compared to the experiments and atomistic simulations in order to validate the CG model ( Table 1). The system density was found to be in good agreement, with the values being about 19% and 27% higher than in the experiments [80] and previous united-atom MD (UAMD) simulations [38], respectively. The thermal expansion coefficient was found to be 52% and 19% higher than in the experiments [81,82] and UAMD [38], respectively. In addition, the squared end-to-end distance (R 0 ), the squared radius of gyration (R g ), and polymer expansion factor (< R 2 0 >/< R 2 g >) of the cis-PI in the melts are within 5% from the UAMD data [38]. R 0 is defined as the average distance between the first and the last CG bead of the cis-PI chain. R g was computed using.
where r i is the position of particle i and r cm is the position of the center of mass of the cis-PI chain. However, the bulk modulus was 3.8-fold lower than in the experiments [83] and 2.6-fold below the previous UAMD simulations [38]. One of the major contributors to this decrease of bulk modulus is the reduction in the degrees of freedom upon coarse-graining monomers into superatom beads, which diminishes the friction between branching atoms and changes entropy [84,85]. Higher mobility of the CG chains causes large fluctuations of the simulation box volume compared to the UAMD simulations [38]. Table 1. Macroscopic and structural properties of cis-PI in melts. R 0 , R g , and < R 2 0 >/< R 2 g > are the end-to-end distance, the radius of gyration, and polymer expansion factor, respectively.

Properties
Exp.

Effect of C 60 Concentration on Macroscopic and Structural Properties of cis-PI Composites
Adding C 60 into the cis-PI composites leads to a linear increase in density ( Figure 1a) due to the mass of the C 60 filler. Figure 1b shows the enhancement of the bulk modulus as a function of the C 60 concentration. A slow increase was found for [C 60 ] < 16 phr (2.17 × 10 −4 GPa per 1 phr), followed by a rapid increase at [C 60 ] = 32 phr (4.81 × 10 −4 GPa per 1 phr). The increases in the density and bulk moduli as a function of the C 60 concentration is similar to previous UA simulations [86]. In addition, the thermal expansion coefficient (γ) and specific heat capacity (c p ) at 300 K were computed using Equations (1) and (2). In the composites, both quantities increased when the amount of C 60 increased (Figure 1c,d), similar to what has been observed in previous MD [87] and experimental studies [88] of polymer carbon-based nanocomposites.  (3)), (c) thermal expansion (Equation (1)), and (d) heat capacity (Equation (2)) of cis-PI and C60 composites as a function of C60 concentration.

Effect of C60 Concentration on Microscopic Properties of cis-PI Composites
To investigate the interaction free energy of C60-C60 and C60-cis-PI in the composites at different C60 concentrations, the potential of mean force (PMF) was calculated from the radial distribution function (RDF; g(r)) [90,91] as:  (3)), (c) thermal expansion (Equation (1)), and (d) heat capacity (Equation (2)) of cis-PI and C 60 composites as a function of C 60 concentration.
The structural changes of the cis-PI chains at different C 60 concentrations were also investigated by monitoring the R 0 and R g . The averages of the R 0 and R g for the last 5 µs  Table 2. For the cis-PI chains, they increased slightly upon increasing the amount of C 60 . The polymer chain expansion factor (< R 0 2 >/< R g 2 >) of the cis-PI in the melts was 6.11, which is in good agreement with previous UAMD simulations (6.20 [38] and 5.44 [77]) and refers to a random coil structure of a cis-PI chain [89]. Table 2. Average end-to-end distance (R 0 ), radius of gyration (R g ) of cis-PI chains, and their polymer chain expansion factor (< R 2 0 >/< R 2 g >) as a function of C 60 concentration.

C 60 Concentration
(phr) To investigate the interaction free energy of C 60 -C 60 and C 60 -cis-PI in the composites at different C 60 concentrations, the potential of mean force (PMF) was calculated from the radial distribution function (RDF; g(r)) [90,91] as: where k B is the Boltzmann constant and T is the absolute temperature. The dimerization free energies of C 60 in cis-PI composites at different C 60 concentrations are shown in Figure 2a. The first, second, and third minima of the C 60 -C 60 PMF corresponding to the distances between the neighboring C 60 molecules are at 1.0 nm, 1.5 nm, and 1.9 nm, respectively. The second minimum is at the lowest free energy, suggesting energetic preference of C 60 dimerization with the cis-PI chain insertion. That the polymer chains fill the space between the C 60 molecules has also been observed in previous studies [38,59,60,92]. The PMF for C 60 and cis-PI interactions is shown in Figure 2b. This result demonstrates the C 60 -cis-PI composites' preference to a dispersed structure corresponding to the snapshots in (Figure 3), which show the dispersed structure for all the C 60 concentrations.

Diffusion of C60 and cis-PI in the Composites
The diffusion coefficients of the C60 and cis-PI in the melts and in the composites w analyzed at different C60 concentrations using the mean squared displacement (MSD) where D is the diffusion coefficient and t is time.
In agreement with previous studies [38,87], the presence of C60 caused a drama decrease in both the C60 and cis-PI diffusion coefficients ( Figure 4). The slow movemen high C60 concentration led to confinement, with the C60 molecules being more mobile th the cis-PI ones. The difference between the diffusion coefficients of C60 and cis-PI creased upon increasing the C60 concentration; the increased interactions between C60 a cis-PI cause a reduction in the diffusive motions of C60 and cis-PI, which is the reason w the bulk modulus of the composites at high C60 concentrations becomes enhanced. N that the diffusion coefficients, extracted from the CG simulations, were one order of m nitude faster than those obtained from the UA simulations [38]. This is due to the red tion of degrees of freedom in the CG model [93]. In Figure 2a,b, the first and second local minima are at 0.8 nm and 1.2 nm, respectively. Regarding the first local minima in Figure 2a,b, the free energy for C 60 -cis-PI is lower than for C 60 -C 60 , indicating a preference for the C 60 -cis-PI interaction. Moreover, in the case of C 60 -C 60 interaction, a concentration dependence was observed in which the free energy decreased slowly upon increasing the C 60 concentration from 4 to 8 phr. For the C 60 -cis-PI, only a minor concentration dependence is present.

Diffusion of C 60 and cis-PI in the Composites
The diffusion coefficients of the C 60 and cis-PI in the melts and in the composites were analyzed at different C 60 concentrations using the mean squared displacement (MSD) where D is the diffusion coefficient and t is time.
In agreement with previous studies [38,87], the presence of C 60 caused a dramatic decrease in both the C 60 and cis-PI diffusion coefficients (Figure 4). The slow movement at high C 60 concentration led to confinement, with the C 60 molecules being more mobile than the cis-PI ones. The difference between the diffusion coefficients of C 60 and cis-PI decreased upon increasing the C 60 concentration; the increased interactions between C 60 and cis-PI cause a reduction in the diffusive motions of C 60 and cis-PI, which is the reason why the bulk modulus of the composites at high C 60 concentrations becomes enhanced. Note that the diffusion coefficients, extracted from the CG simulations, were one order of magnitude faster than those obtained from the UA simulations [38]. This is due to the reduction of degrees of freedom in the CG model [93].

Effect of C60 Concentration on C60 Solvation Free Energy
The solvation free energies of C60 in water, cis-PI melt, and cis-PI-C60 matrix composites were estimated by thermodynamic integration (TI) [78]. The solvation free energy of C60 in water was −99.53 ± 0.03 kJ/mol, which is in good agreement with our previous CGMD calculation using umbrella sampling (−92.6 kJ/mol) [94]. Figure 5 shows the concentration dependence of the solvation free energy. The lowest value (~−115.53 kJ/mol) was found when an added C60 was solvated in a cis-PI melt. It indicates that the cis-PI melt is more favorable than water. Increasing the amount of C60 concentration led to lesser solvation in the composites. Our results suggest that the C60 prefers to interact with the cis-PI chains rather than the other C60 molecules in the cis-PI-C60 composites (Figure 2a,b).

Effect of C60 Concentration on Glass Transition Temperature (Tg)
As discussed above, the CGMD simulations of cis-PI-C60 composites using our model are able to reproduce the results from previous experiments and atomistic MD simulations [50,95]. The glass transition temperature (Tg) is one of the most important parameters determining the physical state and final mechanical properties of polymer composites

Effect of C 60 Concentration on C 60 Solvation Free Energy
The solvation free energies of C 60 in water, cis-PI melt, and cis-PI-C 60 matrix composites were estimated by thermodynamic integration (TI) [78]. The solvation free energy of C 60 in water was −99.53 ± 0.03 kJ/mol, which is in good agreement with our previous CGMD calculation using umbrella sampling (−92.6 kJ/mol) [94]. Figure 5 shows the concentration dependence of the solvation free energy. The lowest value (~−115.53 kJ/mol) was found when an added C 60 was solvated in a cis-PI melt. It indicates that the cis-PI melt is more favorable than water. Increasing the amount of C 60 concentration led to lesser solvation in the composites. Our results suggest that the C 60 prefers to interact with the cis-PI chains rather than the other C 60 molecules in the cis-PI-C 60 composites (Figure 2a,b).

Effect of C60 Concentration on C60 Solvation Free Energy
The solvation free energies of C60 in water, cis-PI melt, and cis-PI-C60 matrix composites were estimated by thermodynamic integration (TI) [78]. The solvation free energy of C60 in water was −99.53 ± 0.03 kJ/mol, which is in good agreement with our previous CGMD calculation using umbrella sampling (−92.6 kJ/mol) [94]. Figure 5 shows the concentration dependence of the solvation free energy. The lowest value (~−115.53 kJ/mol) was found when an added C60 was solvated in a cis-PI melt. It indicates that the cis-PI melt is more favorable than water. Increasing the amount of C60 concentration led to lesser solvation in the composites. Our results suggest that the C60 prefers to interact with the cis-PI chains rather than the other C60 molecules in the cis-PI-C60 composites (Figure 2a,b).

Effect of C60 Concentration on Glass Transition Temperature (Tg)
As discussed above, the CGMD simulations of cis-PI-C60 composites using our model are able to reproduce the results from previous experiments and atomistic MD simulations [50,95]. The glass transition temperature (Tg) is one of the most important parameters determining the physical state and final mechanical properties of polymer composites

Effect of C 60 Concentration on Glass Transition Temperature (T g )
As discussed above, the CGMD simulations of cis-PI-C 60 composites using our model are able to reproduce the results from previous experiments and atomistic MD simulations [50,95]. The glass transition temperature (T g ) is one of the most important parameters determining the physical state and final mechanical properties of polymer composites [96]. To gain more understanding on the influence of the filler concentration on composites' properties, T g was estimated using the change of slope in density (ρ) versus temperature (T) curve [75][76][77]. The ρ-T curves of the cis-PI in the melts and composites are shown in Figure S3. The calculated T g was 201 K, in good agreement with the all-atom (209 K) [77] and united-atom models (223 K) of Sharma et al., [77], and experimental data (200 K) of Brandrup et al. [82] Finally, the concentration dependence of the T g in the composites was determined. The data in Figure 6 show a linear increase with a slope of 0.63 K/phr. This same trend has been reported in previous experiments [34] and simulations [41,87].
Polymers 2021, 13, x 9 of 13 [96]. To gain more understanding on the influence of the filler concentration on composites' properties, Tg was estimated using the change of slope in density (ρ) versus temperature (T) curve [75][76][77]. The ρ-T curves of the cis-PI in the melts and composites are shown in Figure S3. The calculated Tg was 201 K, in good agreement with the all-atom (209 K) [77] and united-atom models (223 K) of Sharma et al., [77], and experimental data (200 K) of Brandrup et al. [82] Finally, the concentration dependence of the Tg in the composites was determined. The data in Figure 6 show a linear increase with a slope of 0.63 K/phr. This same trend has been reported in previous experiments [34] and simulations [41,87].

Conclusions
We have developed a MARTINI force field version 2.1-based [52] parameterization and performed systematic studies of cis-PI in melts and cis-PI-C60 composites using it. We validated the thermodynamic, macroscopic, and microscopic properties of the cis-PI in the melts, and the results show good agreement with prior experimental [80][81][82][83] and computational studies [38,77,94]. Perhaps most surprisingly, only a 0.5% difference with respect to the experimentally found glass transition temperature (Tg) [82] in the cis-PI in the melts was found.
The properties of the NR-fullerene composites at different C60 concentrations (0, 4, 8, 16, and 32 phr) were studied by using CGMD simulations over 200 microseconds. The density, bulk modulus, thermal expansion, heat capacity, and glass transition temperature increased upon increasing the C60 concentration. A slight increase in the R0 and Rg of the cis-PI chains was found for the NR-C60 composites. The interaction energies between the C60 and cis-PI decreased when the amount of C60 increased. The diffusion of the cis-PI and C60 was slowed down by increasing the C60 concentration. The decrease of the interaction free energies and the diffusion coefficients resulted in an enhancement of the bulk modulus of the composites. Moreover, the solvation free energies of the C60 in the cis-PI matrix composites increased upon increasing the C60 concentration. These results suggest that C60-cis-PI interactions are more preferable than C60-C60 self-interactions.
The CG model introduced here allows simulations of large systems over long periods of time to study properties such as the impact of the filler concentration on rubber composites at the molecular level. This validated cis-PI-C60 CG model based on the MARTINI force field will also enable simulations of the advanced rubber materials and their interactions with biological molecules.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1: Details of the CG force field for cis-PI, Table S2: Number of molecules and simulation time, Table   Figure 6. Glass transition temperature (T g ) of cis-PI-C 60 composites at different C 60 concentrations.

Conclusions
We have developed a MARTINI force field version 2.1-based [52] parameterization and performed systematic studies of cis-PI in melts and cis-PI-C 60 composites using it. We validated the thermodynamic, macroscopic, and microscopic properties of the cis-PI in the melts, and the results show good agreement with prior experimental [80][81][82][83] and computational studies [38,77,94]. Perhaps most surprisingly, only a 0.5% difference with respect to the experimentally found glass transition temperature (T g ) [82] in the cis-PI in the melts was found.
The properties of the NR-fullerene composites at different C 60 concentrations (0, 4, 8, 16, and 32 phr) were studied by using CGMD simulations over 200 microseconds. The density, bulk modulus, thermal expansion, heat capacity, and glass transition temperature increased upon increasing the C 60 concentration. A slight increase in the R 0 and R g of the cis-PI chains was found for the NR-C 60 composites. The interaction energies between the C 60 and cis-PI decreased when the amount of C 60 increased. The diffusion of the cis-PI and C 60 was slowed down by increasing the C 60 concentration. The decrease of the interaction free energies and the diffusion coefficients resulted in an enhancement of the bulk modulus of the composites. Moreover, the solvation free energies of the C 60 in the cis-PI matrix composites increased upon increasing the C 60 concentration. These results suggest that C 60 -cis-PI interactions are more preferable than C 60 -C 60 self-interactions.
The CG model introduced here allows simulations of large systems over long periods of time to study properties such as the impact of the filler concentration on rubber composites at the molecular level. This validated cis-PI-C 60 CG model based on the MARTINI force field will also enable simulations of the advanced rubber materials and their interactions with biological molecules.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/polym13224044/s1, Table S1: Details of the CG force field for cis-PI, Table S2: Number  of molecules and simulation time, Table S3: Solvation free energy of an added fullerene in water, cis-PI in melt and composite at different C60 concentrations, Figure S1: Time evolution of end-to-end distance (R0) with varies C60 concentration from 0 to 32 phr, Figure S2: Time evolution of radius of gyration (Rg) with varies C60 concentration from 0 to 32 phr, Figure S3: Density versus temperature of cis-PI-C60 composites.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.