Molecular Dynamics Simulations of Fullerene Diffusion in Polymer

Article pubs.acs.org/Macromolecules

Molecular Dynamics Simulations of Fullerene Diffusion in Polymer Melts Igor V. Volgin,† Sergey V. Larin,† Enrique Abad,‡ and Sergey V. Lyulin*,†,§ †

Institute of Macromolecular Compounds, Russian Academy of Sciences, Bolshoj pr. V.O., 31, 199004 Saint Petersburg, Russia Departamento de Física Aplicada and Instituto de Computación Científica Avanzada (ICCAEX), Centro Universitario de Mérida, Universidad de Extremadura, E-06800 Mérida, Spain § Physical Faculty, Saint-Petersburg University, Ulyanovskaya str. 1, 198504 Petrodvorets, Russia ‡

ABSTRACT: We employ all-atom molecular dynamics simulations up to microsecond time scales to study diffusion of fullerene nanoparticles (C60 and its derivative PCBM) in a polyimide matrix above its glass transition temperature. A detailed examination of the fullerene mobility in the embedding polymer system reveals the presence of different diffusion regimes (ballistic, subdiffusive, and normal diffusive). The microscopic origin of the observed subdiffusive regime is discussed by comparing the behavior to the one displayed by different anomalous diffusion processes, namely, continuous time random walk (CTRW), random walk on a fractal (RWF), and fractional Langevin equation (FLE). A series of statistical tests suggests that the FLE framework is the more appropriate one to describe subdiffusion of fullerenes in our system. Furthermore, a comprehensive analysis of the self-part of the van Hove function shows that the normal diffusion regime observed at long times displays a nonclassical behavior characterized by the simultaneous presence of several Gaussian peaks. We ascribe this behavior to a mechanism of diffusion by hopping. Until recently, it was commonly believed that hopping of tracer particles in polymer systems is only relevant for particle sizes that are of the order of the reptation tube diameter (dT), while smaller particles are considered to slip through the entanglement mesh. However, our results provide direct evidence for hopping as a relevant mechanism for diffusion of particles whose sizes are commensurate with the correlation length (ξ) of the polymer system. These results emphasize the importance of local interactions at the atomic level for the understanding of nanoparticle dynamics in polymer melts.

I. INTRODUCTION Incorporating nanoparticles into polymer matrices as fillers is a common way to produce polymer composite materials (PCMs). Aromatic polyimides (PI) constitute a promising class of polymers which can be used as matrices in highperformance PCMs. This is due to their unique characteristics, such as thermal stability, thermoplasticity, and chemical resistance, all of which are of crucial importance for practical applications of PI-based composites in modern industry.1,2 Different types of nanoparticles can be used as fillers in PIbased nanocomposites, thus allowing one to obtain materials with improved characteristics.3−8 However, to achieve a better performance of such materials during its processing by melt extrusion, the filler species needs to be dispersed in a rather uniform way throughout the polymer matrix. In general, this is difficult to achieve because of the ongoing nanoparticle aggregation processes.9,10 In this context, the improvement of the currently available processing techniques with the aim of optimizing the molecular design of advanced nanocomposite materials requires a deep understanding of the dynamical behavior of nanoparticles in polymer melts in general and of the underlying mechanisms of diffusive transport in particular. © XXXX American Chemical Society

Generally speaking, the diffusion of nanoparticles in polymer systems is controlled by different factors such as hydrodynamic forces, the segmental motion of the polymer chains, and topological constraints.11−17 The dynamic properties of the filler particles in nanocomposites are considerably influenced by several key factors, namely, nanoparticle structure (e.g., shape and type of surface modification), molecular weight and topological architecture of the polymer chains, nature of the polymer−particle interactions, and ratio between the characteristic size of the particle and that of the polymer chains. The latter factor plays a decisive role in the study of nanoparticle diffusion by means of the existing experimental methods.18−31 As a matter of fact, most of the aforementioned works are devoted to nanoparticle diffusion in polymer solutions,18−27 whereas those dealing with nanoparticle mobility in polymer melts remain scarce.28−31 These studies provide a direct evidence for deviations of the nanoparticle mobility from the Stokes−Einstein relation, thus confirming the presence of the Received: September 20, 2016 Revised: February 6, 2017

A

DOI: 10.1021/acs.macromol.6b02050 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules scale-dependent viscosity effect32 predicted and described by various theoretical frameworks such as mode-coupling theory,11 statistical dynamics approaches,12−14 or scaling theory.15−17 In particular, according to the scaling theory recently developed by Cai et al.,17 there are three different scenarios of particle motion in the polymer matrix depending on the ratio between the particle size and the main structural length scales in the polymer system, such as correlation length (ξ) or reptation tube diameter (dT). As it turns out, the dynamics of small particles (R < ξ) does not depend on the structural relaxation of the polymer medium, while for intermediate particle sizes (ξ < R < dT) the diffusion process begins to be affected by the segmental motion of the polymer chains. Finally, diffusion of comparatively large particles (R > dT) is controlled to a large extent by the relaxation dynamics of the entanglement mesh. By measuring the diffusion coefficient of the particles and the viscosity of the polymer system, one finds that the polymer bulk viscosity only has a significant influence on the diffusion of sufficiently large particles, whereas particles with a diameter smaller than dT tend to diffuse faster than predicted by the Stokes−Einstein relation. A possible explanation is that a mechanism of so-called “hopping diffusion”14,17 may be at play, an issue which was quantitatively addressed in the framework of various statistical dynamics theories.13 According to the theoretical predictions obtained so far, hopping is relevant only for particle sizes 2R/dT ∼ 1.5−2. To date, even though much effort has gone into testing these predictions,20,24,29 the direct observation of the dynamic properties of medium-sized particles in experiments and simulations still remains challenging, especially for particles whose size is commensurate with the typical correlation length of polymer melts. On the other hand, the diffusion properties of such particles are of great importance for the molecular design of novel PCMs having them as fillers. In the above context, computer simulations are a powerful tool to study the diffusion of medium-sized particles. Having said this, most simulations concerning polymer−nanoparticle systems were performed in the framework of coarse-grained or bead−spring models33−49 and aimed at testing fundamental ideas of nanoparticle or polymer diffusion. These works provide clear evidence for the fact that nanoparticle dynamics is controlled by a ratio of the structural length scales in the system,33−36 but it was also found that overall nanocomposites properties depend significantly on the particular features of the “polymer−nanoparticle” interactions45−49a distinctive feature of each nanocomposite system. Therefore, dynamic properties of particles in polymer systems should be investigated by means of simulations that take into account all the structural details of the nanocomposite system components, along with a correct parametrization of the “polymer−nanoparticle” interactions. Until recently, diffusion in PI-based systems had only been studied in PI gas separation membranes,50−59 where gas diffusion is determined mainly by free volume fluctuations of the system. However, polyimides may also be considered as interesting candidates to investigate the diffusive transport of embedded medium-sized particles, e.g., fullerenes. For these systems, the main role in particle diffusion is played by the nanoviscosity of the polymer, which in turn may lead to interesting effects. Nevertheless, there is only a limited number of works on fullerene-containing systems, and their main focus is on elucidating the effect of fullerenes either on the local segmental mobility of the polymer chains or on diffusion of small gases in nanocomposites.60,61 In contrast, the dynamic

properties of the fullerenes themselves have only been investigated in biological systems.62−66 However, diffusion of fullerenes in thermoplastic polymers is of great interest, since the introduction of fullerenes may increase thermal stability due to radical scavenging processes,67 and are known to influence wear resistance68 as well as optical69 and dielectric70 characteristics of PCMs. Because of such effects, fullerene-containing PCMs have proven to be useful for optical systems, nanoelectronics devices, and photovoltaic applications.8 For the present study, we performed all-atom molecular dynamics (MD) simulations of fullerene nanoparticles to mimic their diffusion up to microsecond time scales in PI matrices above the glass transition temperature. Our goals are to investigate the influence of the chemical structure of fullerene on its mobility in the PI melt as well as to test a number of theoretical considerations on the underlying mechanisms of nanoparticle diffusion in the framework of models incorporating atomistic details. The outcome of such studies may provide new opportunities to control nanoparticle dispersion and aggregation in nanocomposite systems as well as stimulate the development of new methods for the control of morphology in nanocomposite materials.

II. MATERIALS AND METHODS We studied key diffusion properties of the fullerene C60 and its derivative [6,6]-phenyl-C61-butyric acid methyl ester (PCBM) in the matrix of thermoplastic PI R-BAPB, synthesized from 1,3-bis(3′,4dicarboxyphenoxy)benzene (dianhydride R) and 4,4′-bis(4″-aminophenoxy)diphenyl (diamine BAPB), and characterized by a comparatively high glass transition temperature Tg = 479 K71 (Figure 1).

Figure 1. Chemical structure of the R-BAPB polyimide, C60 fullerene, and PCBM. A segment comparable in size with the fullerene diameter is marked by the red arrow (see text for details). Computer simulations were carried out by means of the MD method implemented in the Gromacs 4.6 package.72,73 To perform our MD simulations, we used a computational scheme together with simulation parameters that had already been successfully applied in our previous research74−81 on physical properties of bulk polyimides,74−79 PI-based nanocomposites containing bare80 and modified nanotubes,81 and graphene.82,83 In the present work, the cubic cell with periodic boundary conditions contained 27 PI chains with a degree of polymerization (DP) equal to 8 as well as a single nanoparticle (either C60 or PCBM). The characteristic sizes of the simulated systems are given in the Table 1.

Table 1. Characteristic Sizes of the Simulated Systems R-BAPB polyimide

B

box size (nm)

fullerene diameter (nm)

Rg (nm)

lp (nm)

monomer unit (nm)

entanglement length (nm)

packing length (nm)

∼6

∼1

∼3.7

∼1.3

3.4

∼10

∼0.1


Article

Macromolecules

details of their chemical structure, we first evaluated the meansquared displacement (MSD) of the fullerene center-of-mass:

The main reason for studying systems with DP = 8 was to eliminate the possible influence of the entanglement of polymer chains on the final results. Thus, we used DP = 8, corresponding to the onset of the polymer regime (approximately three entanglements per chain). Additional simulations of systems consisting of longer polymer chains (DP = 16) were also performed to study the influence of the DP value on the obtained results. As expected, upon such an increase of the polymerization degree the simulation times needed to obtain equilibrated polymer samples become much longer. However, the results for fullerene diffusion in this system during equilibration are almost the same as in composites based on R-BAPB polyimide with DP = 8. This result suggests that further increasing of the polymerization degree of the system does not lead to any significant differences in nanoparticle diffusion and thus has no effects on the main results of the present work. The simulations were performed in the NpT ensemble at a pressure of 1 bar and a temperature of 600 K, which lies in the range between the glass transition temperature of the PI bulk samples and the threshold value for its thermal destruction.74 Atomic partial charges were taken to be zero, as their values are relatively small compared to the significant electrostatic interactions which may be induced by the highly polarized groups contained by the polyimide structure (for example, sulfone).74 On the other hand, since fullerenes are charge neutral molecules, simulations without partial charges will strongly decrease the time cost of the runs without loss of the specific “particle−polymer” interactions. Moreover, it was recently shown that under the aforementioned simulation conditions no viscoelastic effects are observed in the polyimide matrix, implying that it can be considered a pure polymer melt.78 Finally, in order to examine the influence of the fullerene addition on the equilibration time of the system, we estimated the time evolution of the radius of gyration Rg of R-BAPB PI chains (Figure 2).

⟨Δr 2(Δt )⟩ = ⟨(r(t0 + Δt ) − r(t0))2 ⟩t0

(1)

where ⟨...⟩ denotes a moving time average (the subindex t0 indicates explicitly the parameter with respect to which the average is taken). The resulting curves (Figure 3) are characterized by the presence of several diffusion regimes which can be

Figure 3. A log−log plot of the center-of-mass MSD of the fullerenes in the R-BAPB polyimide melt. The inset shows the MSD in linear coordinates. Dashed lines indicate the slopes predicted by scaling theory for medium-sized particles (ξ < R < dT).

distinguished by the different slopes of the MSD in logarithmic coordinates: ballistic (MSD ∼ Δt2, Δt < 0.5 ps), subdiffusive (MSD ∼ Δt0.4, 10−3 < Δt < 10 ns), and normal diffusion (MSD ∼ Δt, Δt > 20 ns) regimes. In the ballistic and subdiffusive regimes the slopes of the curves (as well as their typical spatial and temporal scales) are independent of the surface modification of the fullerenes. This fact can be accounted for by the insignificant influence of the grafted “tail” in the PCBM on its effective size. While the ballistic regime and the normal diffusion regime are both well-known from the theory of Brownian diffusion, the onset of an intermediate subdiffusive regime is challenging and calls for an extension of the standard theory. In the theory of anomalous diffusion processes, one speaks of “subdiffusion” when the MSD displays a sublinear time dependence (MSD ∼ Δta, 0 < a < 1) in the long time limit. In practice, this applies to time scales significantly larger than the duration of the experiment or significantly larger than the typical relaxation times of the system at hand. However, in our case subdiffusion does not persist during the entire duration of the numerical experiment; hence, the term “transient subdiffusion”a typical feature of polymer systemsis appropriate here.85 In particular, according to predictions based on scaling theory,15−17 subdiffusion may be brought about by relaxation processes of polymer chain segments commensurate with the size of the diffusing particles. In order to test this theoretical assumption for our systems, we compared the MSD of the C60 fullerene with the MSD of commensurate R-BAPB segments (Figure 4). This PI fragment was chosen because of its nearly flat configuration, which reduces significantly the impact of the structural fluctuations on the calculated characteristics. From the figure one can see the coincidence between the MSDs of

Figure 2. Time dependence of the radius of gyration Rg corresponding to R-BAPB polyimide chains for systems containing fullerene C60 (black circles) and PCBM (red squares) as well as for the bulk system (orange triangles). The theoretical value of Rg is given by the dashed line.

As can be seen from the figure, Rg approached its upper limit after approximately 1 μs, as in the case of the unfilled samples.74−79 Thus, the addition of the fullerenes did not have a significant impact on the evolution of the Rg values. The obtained limiting values of the PI RBAPB radius of gyration are in good agreement with theoretical calculations based on the virtual-bond formalism,84 which yield Rg = 3.7 nm.74 Hence, our subsequent studies of the nanoparticle dynamics were carried out on fullerene trajectories obtained after the characteristic equilibration time of 1 μs for each system.

III. RESULTS AND DISCUSSION A. Mean-Squared Displacement. In order to investigate the dependence of dynamic properties of fullerenes upon the C


Article

Macromolecules

the guidelines given by Meroz and Sokolov.88 To begin with, we tested the stationarity of the displacement increments as well as the ergodicity of the underlying process. Stationarity was tested by means of the quantile line test88 performed on an ensemble of five 100 ns trajectories. These trajectories were obtained by carrying out additional simulations from different points of the main microsecond trajectory. The obtained quantile lines, each of them calculated from its own trajectory, are depicted in Figure 5. The obtained

Figure 4. A log−log plot of the center-of-mass MSDs of the C60 fullerene (black line), the R-BAPB chain (blue), and an R-BAPB segment (red) comparable in size with the fullerene (Figure 1).

the fullerene and the chosen segments on intermediate times scales of 10−3−10 ns corresponding to the subdiffusive time window. At Δt ∼ 10 ns both curves separate, precisely at the end of the subdiffusion regime. This behavior suggests that for intermediate time scales the fullerene motion is coupled to the translational motion of the R-BAPB segments. Uncoupling appears to take place at time scales of ca. 10 ns, which corresponds to a fullerene displacement of about 1 nm, a quantity of the order of the nanoparticle diameter and coincident with the polymer persistence length. This effect is observed for both fullerenes; significant differences between their MSDs only manifest themselves on sufficiently long time scales, especially in the normal diffusion regime. One can see that the grafted “tail” slows down PCBM diffusion with respect to C60 diffusion. This may be due to the presence of additional geometrical restrictions hindering fast diffusive displacement. In view of the qualitative agreement of the MSD curves for PCBM and C60 in the ballistic and subdiffusive regimes, we shall restrict our subsequent discussion of the observed subdiffusive regime to the C60 fullerene so as to avoid redundancies of the available data. B. Theoretical Framework for the Subdiffusive Regime. In contrast to scaling theory, our results show that the slope of the MSD in the subdiffusive regime is 0.02 ns) we compared the MSD obtained from five 200 ns trajectories with the MSD from a single 1000 ns trajectory. As a result, we found that the MSD curves overlap, implying that the trajectory-to-trajectory scattering is negligible and that our diffusion process is ergodic. D


Article

Macromolecules To confirm the conclusions provided by the stationarity test and the ergodicity test, we carried out an additional test whose aim is to distinguish nonergodic CTRW processes from ergodic processes in the FLE class, namely, the p-variation test.89 The main idea is to divide the trajectory into 2N intervals (tj, tj+1) [where N is a parameter of calculations] and to plot partial sums V(p) N given by the expression 2N − 1

V N(p)(t )

=

∑

|r(min{t j + 1 , t }) − r(min{t j , t })| p (2)

j=0

V(p) N ’s (p)

against the time t until which the sum runs, where the are the Nth-order approximants of the p-variation V (t) = limN→∞V(p) N (t). According to the study of Magdziarz et al.,89 for an FBM process V(2) N tends to infinity with increasing N, while for a CTRW process V(2) N should stabilize. In our calculations we analyzed one 50 ns trajectory recorded in 10−5 ns steps. By increasing value of N from 4 to 17, we found that the slopes of the corresponding V(2) N curves increase (Figure 7), suggesting

Figure 8. (a) Non-Gaussian parameter α2 as a function of time. Dashed lines delimit the subdiffusion interval. (b) Corresponding values of the MME ratio.

trajectory. In the subdiffusive interval there are some deviations from the Gaussian form. However, direct fitting of the DD with Gaussian curves (Figure 9) suggests that such deviations are not very significant.

Figure 7. VN(2) curves as a function of N obtained from a single 50 ns trajectory. Figure 9. Displacement distribution (DD) functions 4πr2GS(r,Δt) for various subdiffusive time intervals Δt (symbols). Solid lines represent corresponding Gaussian fits.

that V(2) N tends to infinity with increasing N. The above results exclude a CTRW process as the underlying subdiffusive process, thereby confirming our previous conclusion based on the outcome of the stationary test and the ergodicity test. Our next step was to test the gaussianity of the distribution of displacements (DD) in the time interval corresponding to the subdiffusive regime. This distribution can be easily obtained from the self-part of the van Hove function (VHF) Gs(r , Δt ) = ⟨δ[r − (|r(t0 + Δt ) − r(t0)|)]⟩t0

To support this idea, we recall that in addition to the nonGaussian parameter some authors100 use the minimal excursion ratio (MME) to discriminate between different subdiffusive models:

(3)

MME =

by multiplying it by a factor 4πr2. As a result, the DD function (or 4πr2Gs(r,Δt)) gives the probability density to find a particle at a distance r from the initial position after a time interval Δt. Figure 8a shows the non-Gaussian parameter α2, defined by the following equation: α2(Δt ) =

3⟨Δr 4(Δt )⟩ −1 5⟨Δr 2(Δt )⟩2

⟨Δrmax 4(Δt )⟩ ⟨Δrmax 2(Δt )⟩2

(5)

In our case the regular MME ratio is about 1 (see Figure 8b), thereby suggesting that the subdiffusive transport process is compatible with FBM theory.88 Inasmuch as the DD functions are not perfectly Gaussian, it is reasonable to ask whether a RWF process could also be a suitable candidate to describe our simulation results, especially in view of the fact that the analytical form of the DD function for a RWF process is not known in general.88 To discern whether an RWF or an FLE/FBM process underlies the subdiffusion process, we tested how the fullerene trajectory fills the space by means of the so-called growing sphere test.91 The main idea behind this test is to calculate the probability for a particle to remain inside a sphere of radius r0 whose growth rate

(4)

This parameter quantifies the degree of Gaussianity of a given probability distribution (a zero value provides strong evidence for a perfect Gaussian distribution). For short time intervals (Δt < 0.2 ns) α2 was obtained by ensemble averaging over ten 50 ns trajectories, while for long times (Δt > 0.2 ns) it was evaluated from a single 1000 ns E


Article

Macromolecules is determined by the slope a of the MSD curve in the subdiffusive regime: Pr(r ≤ r0t a /2) ≈

1 N

N

∑ θ(ri(t ) − r0t a /2) i=1

(6)

where N is the number of trajectories, ri is the distance from the starting point, and θ(r) stands for the Heaviside step function. According to Meroz and Sokolov,88 a suitable value of r0 is given by the average length of a single trajectory step. In our calculations we took r0 = 0.17 nm and a = 0.4 nm2/ns. Besides, 50 subsets of the main microsecond trajectory were chosen, where the length of each subset was equal to 10 ns. The trajectory subsets were nonoverlapping and separated by 10 ns time intervals. The trajectory length of 10 ns was chosen to guarantee that the fullerene remained inside the subdiffusive time window, where its motion is hindered by the cage effect induced by the polymer chain segments (see below). The results provided by the growing-sphere test are displayed in Figure 10. In spite of the fact that the results are affected by

Figure 11. A log−log plot of the absolute value of the displacement correlation (DC) function as a function of the time interval Δt for the C60 fullerene motion in the R-BAPB melt. The inset shows timedependence of the DC function in linear coordinates.

increased (decreased) probability for two consecutive displacements to take place in the same direction. Plotting the modulus of the DC function allows one to identify time scales when the correlation tendencies are changed (Figure 11). The first “peak” (Δt ∼ 10−3 ns) refers to the typical time when the ballistic regime ends. Beyond this time, i.e., on intermediate time scales, the DC function is seen to increase linearly with a slope of 0.4 nm2/ns. Since DC ∼ Δt2−2H, the slope of the DC curve yields the Hurst exponent H of the corresponding FLE. In our case, the value of the slope corresponds to a Hurst exponent of 0.8, which indicates a tendency for the displacements to be anticorrelated in the subdiffusive regime.92 This suggests that at intermediate times the fullerene’s motion is hindered by the cage formed by the chain segments. We close the present subsection by noting that the FLE is usually employed to describe viscoelastic systems, where memory effects are incorporated via power-law functions,93 Mittag−Leffler functions,94 Prony series,95 etc. All the memory kernels are defined via the correlation function of the thermal noise.96 Generally, the type of memory function is deduced from the decay of the velocity autocorrelation function97 via the Mori−Zwanzig approach. Some authors attempted to incorporate memory kernels obtained from the polymer reference interaction site model (PRISM) and the mode coupling theory (MCT) describing the relaxation of polymer medium.11−14 However, in these works only the normal diffusion regime was successfully reproduced, whereas in the present work the FLE provides a reasonably good description of the subdiffusive regime. C. Diffusion in the Long-Time Regime: Onset of NonGaussian Brownian Motion. As already anticipated, it is interesting to note that in the normal diffusive, long time regime the fullerene tends to diffuse faster than the PI chain segments (their MSD curves separate at times around 10 ns; cf. Figure 4). It is also worth mentioning the presence of secondary peaks on the DC curves in the range 10−30 ns (Figure 11), which correspond very accurately to the time scales when the MSDs of the fullerene and the segments intersect, thereby setting an end to the subdiffusive regime. In order to shed light on the microscopic origin of the observed behavior, we analyzed the long-time behavior of the DD function (eq 3). From the theory of conventional Brownian

Figure 10. Probabilities to find the C60 fullerene inside a growing sphere as a function of its initial size r0.

the r0 value, the probability clearly approaches unity for sufficiently long times, implying that the fullerene trajectory fills the available space homogeneously. In view of this result, the RWF process should be excluded as a possible origin of subdiffusion. From the overall picture provided by the different statistical tests, we conclude that the most suitable theoretical framework of fullerene subdiffusion is provided by the FBM model. However, FBM does not account for the ballistic regime of motion. This shortcoming can be corrected by means of the fractional Langevin equation (FLE) which, as already mentioned, is an extended version of the FBM theory.92 In the framework of both theories, the main parameter is the socalled Hurst exponent H, which characterizes the spatial correlations displayed by the transport process as well as the time decay of the memory kernel of the associated FLE. The value of H was computed from the displacement correlation (DC) function92 (Figure 11): DC(Δt ) = ⟨δ r(t0 , Δt )δ r(t0 − Δt , Δt )⟩t0

(7)

where δr(t0,Δt) = r(t0,Δt) − r(t0) is the displacement vector of the fullerene after a time interval Δt. According to its definition, the DC function quantifies the time correlation between consecutive displacements of time duration Δt. Positive (negative) DC values correspond to an F


Article

Macromolecules

identification of fullerene jump events, we first calculated from a 1 μs trajectory the radial distribution function (RDF) referring to the distances between the fullerene center-of-mass and all the polymer atoms in the system. To first approximation, the RDF can be related to the potential of mean force (PMF) describing the interaction between the fullerene and its environment via the formula PMF = −kBT ln(RDF).112 The resulting PMF profile (Figure 13) shows the presence of a potential barrier associated with the first layer of monomers

motion, it is known that in the long-time limit (MSD ∼ Δt) the DD function should follow a Gaussian form with a peak value corresponding to the most probable traveled distance during the time interval Δt. This follows from the gaussianity of the self-part of the van Hove function for Brownian motion:98 GS(r , Δt ) = (4πDΔt )−3/2 e−r

2

/4DΔt

(8)

However, in our case the DD function tends to deviate from the usual unimodal Gaussian form; instead, a multipeaked DD curve (Figure 12) is seen to emerge as Δt increases. Moreover,

Figure 13. Potential of mean force, PMF, along the distance between the fullerene center-of-mass and the polymer atoms in melt. The dotted line indicates the transition distance rmax between two possible fullerene states, i.e., “trapped in the cage” (r < rmax) and “diffused to the bulk” (r > rmax).

Figure 12. DD functions for C60 diffusion in different time scales, i.e., Δt = 10, 30, and 60 corresponding theoretical curves (solid lines). functions are also depicted by dotted lines. Their by the solid black line.

surrounding the fullerene (first coordination shell). This shell forms the cage barrier which the fullerene has to overcome in order to perform a jump. Specifically, the fullerene is considered to have escaped from the cage in case that its displacement is larger than the effective distance rmax = 1.1 nm (Figure 13). This parameter was used to calculate the probability R(t) that a fullerene located inside a monomer cage at t = 0 still remains inside the cage up to time t. One has R (t) = 1 − ⟨p(t)⟩ ≈ e−t/τ, where p(t) = 0 if the displacement r < rmax and p(t) = 1 otherwise. In this equation, τ is the mean residence time of the fullerene in the cage. To perform our analysis, we took the average over a single trajectory of length equal to 1 μs. Strictly speaking, the algorithm that we used does not account for the possibility of reversible escape, i.e., for the fullerene jumping back into the cage. However, these are very rare events and therefore statistically irrelevant in the present context. The resulting correlation function R(t) (Figure 14a) was fitted with an exponential function, from which the value τ ∼ 18 ns was obtained. Remarkably enough, this time matches very precisely the time scale (Δt) corresponding to the decoupling between segmental and fullerene MSDs, which marks the onset of the fullerene hopping regime. In particular, the behavior of the DD functions for time scales from 18 to 30 ns is characterized by the onset of a secondary local peak (or “shoulder”) as Δt increases further. Hence, this additional evidence strongly suggests that the deviations from the Gaussian approximation observed in Figure 14b are likely to be a sign of a first fullerene jump, which it takes after having dwelt inside a cage for approximately 18 ns.

R-BAPB melt for ns (symbols) and Various Gaussian sum is represented

the observed DD curve can be decomposed into several Gaussian functions with different peak positions (Figure 12), as opposed to the DD curves corresponding to the subdiffusive regime (Figure 9). The onset of several peaks suggests that after a given time interval there are several values of the displacement which may occur with a comparatively large probability. Recently, a mechanism of diffusional hopping was invoked to explain a similar effect observed for diffusion in gases and liquid crystals.99−106 Generally speaking, in a situation like the above onewhere the van Hove function (VHF) deviates from the Gaussian form even at times when the MSD remains linearone speaks of “anomalous yet Brownian diffusion”.107 Wang et al. classified the VHF deviations into four families depending on whether its large-r “tail” is exponential or Gaussian and on how fat it is. Various specific interpretations based on different microscopic mechanisms leading to either Gaussian or exponential tails have been put forward, namely, diffusing diffusivity,108 walking confined diffusion,109 spatially heterogeneous dynamics,110 and hopping diffusion.23 Here, we shall adopt the latter one as a starting hypothesis needing additional evidence. To this end, we used a method developed by Laage and Hynes111 to compute the residence time of the fullerene in a transient cage formed by nearby monomers. To establish a clear-cut criterion for the G


Article

Macromolecules

Figure 15. The 100 ns fullerene center-of-mass trajectory of a C60 fullerene in the R-BAPB melt. The arrow indicates a jump event.

scales (Δt < 10−3 ns) the fullerene moves almost freely, and ballistic motion takes place. At longer times, the motion of fullerene is slowed down by geometrical restrictions (cage effect) imposed by chain segments commensurate with the fullerene, resulting in the onset of a subdiffusive regime (10−3 ns < Δt < 10 ns). In this regime, the fullerene follows the motion of the polymer chain segments as the cage formed by the latter progressively relaxes (no plateau region is observed in the corresponding MSD curves). At even longer times, the fullerene eventually performs a long jump. As a result of this hopping event, the fullerene approaches its terminal diffusion regime, corresponding to time scales Δt > 10 ns. As far as the mesoscopic origin of the different regimes is concerned, we showed that a FLE approach provides a unified description for both the ballistic and the subdiffusive regime. In contrast, the FLE does not provide an intuitive description of the hopping events observed in the normal diffusive regime. Hence, additional assumptions are needed to construct an overall theory for the nanoparticle motion in the melt.

Figure 14. (a) Time correlation function R(t) associated with the probability of the fullerene being trapped in the interior of the first coordination shell (r < rmax) of monomers (black line). The red line describes the corresponding exponential fit. (b) Displacement distribution functions for C60 diffusion in R-BAPB melt for time scales Δt = 18, 24, and 30 ns (solid lines) accompanied by the corresponding theoretical Gaussian approximations (dotted lines).

Further, as time goes by the number of local peaks grows along with the number of long jumps made by the fullerene. The distance between adjacent peak locations is nearly the same for all the DD curves and turns out to be of the order of the average jump length distance (∼1.4 nm). Interestingly, the latter is of the same order as the PI persistence length (lp ∼ 1.32 nm for R-BAPB).74 Thus, after the first jump event, the fullerene begins to diffuse faster than the chain segments. We also tracked the fullerene center-of-mass trajectory during a single statistical realization (Figure 15). It shows that fullerene diffusion is characterized by a repeated sequence of short steps localized in a limited region of the space followed by long, almost unidirectional multistep displacements connecting such localized areas. These long displacements may be interpreted as jumps consisting of several successive steps, nearly all of them being aligned along a given direction (Figure 15). Thus, considering the DD curve for Δt = 30 ns, we conclude that short displacements correspond to fullerene diffusion in localized areas, whereas the onset of a second local peak at larger distances is most likely explained by the existence of hopping events. We believe that the physical reason behind the confined motion of the fullerene inside localized regions may be caging caused by segments of the polymer chains, leading to the onset of the subdiffusive regime on sufficiently short time scales (Figure 4). We are now in the position to put forward a plausible picture simultaneously accounting for the different regimes observed for the motion of fullerenes in a polymer melt. On short time

IV. CONCLUSIONS In the present work we employed all-atom MD simulations up to the microsecond time scale to investigate characteristic transport properties of C60 fullerenes as well as of PCBM derivatives diffusing in a polyimide matrix above the glass transition temperature. The present study spans a broad range of length and time scales, thereby comprising different transport regimes. Our computational approach provides a powerful tool to characterize the specifics of the fullerene motion by means of a detailed account of its chemical structure and that of the environment. The analysis of the diffusional behavior of fullerenes up to the microsecond time scale revealed the presence of three different diffusion regimes: ballistic, subdiffusive, and normal diffusive. We found convincing evidence for the hypothesis that the subdiffusive regime is induced by the coupling between the fullerene motion and the translational motion of chain segments commensurate with the size of the fullerene. Even though Cai et al. had already predicted a similar effect for particles smaller than the entanglement mesh size,16 it is H


Article

Macromolecules surprising to find a subdiffusive regime for particles whose size is of the order of the polymer correlation length in the melt. We subsequently applied a number of statistical tests to shed light on the microscopic origin of subdiffusion and to investigate which features of three main model classes CTRW, RWF, and FBMhold. By applying the scheme proposed by Meroz et al.,88 we investigated whether the underlying transport process is stationary and ergodic, whether the associated displacement distribution (DD) function is Gaussian, and whether the fullerene trajectory fills the space homogeneously. In spite of the fact that the DD is not perfectly Gaussian in the subdiffusive regime, our tests indicate that a fractional Langevin equation (FLE) with Hurst parameter H = 0.8 provides a reasonably good description of the fullerene motion. Physically, this value corresponds to the case where consecutive displacements in opposite directions are more probable than consecutive displacements in the same direction. In our systems, this feature arises from the motion of the fullerene in the cage formed by the nearest polymer segments, implying that the cage exerts a viscoelastic effect on the fullerene. For this reason, we believe that the use of the term “nanoviscoelasticity” is appropriate in the present context, in spite of the previously known fact that as shown previously, the polymer melt is not viscoelastic on a macroscopic scale. This observation may provide new opportunities to characterize mechanical properties of the polymer melts at nanoscale, a line of research which we aim to pursue in future works. In the normal diffusive regime, our analysis was again based on the form of the corresponding DD functions. The latter were seen to develop multiple Gaussian peaks in the course of time, thereby supporting the hypothesis that the long-time diffusion of the fullerene is governed by hopping events. Further support for this assumption comes from the calculation of the mean residence time of the fullerene inside the cage as well as from direct observation of the fullerene trajectory. As a result, the average mean residence time was found to match the typical time scale for decoupling between the motion of polymer segments and the fullerene. Moreover, DD functions corresponding to time scales larger than the mean residence time were found to differ from Gaussian form due to the onset of a first local peak, associated with the first jump event. Thus, the distance between the locations of the adjacent peaks displayed by the DD functions was interpreted as an average jump length. Interestingly, the latter was found to be of the same order as the polymer persistence length. Of special interest is the discussion of the above finding from the standpoint of modern theories of nanoparticle diffusion in polymer melts, such as scaling theory17 and statistical dynamical theory.14 According to these theories, hopping is relevant only for particles with sizes 2R: dT < 2R < dc (where dc is some crossover particle size and dT is the reptation tube diameter), corresponding to the regime 2R/dT ∼ 1.5−2. Diffusion of larger particles is dominated by polymer chain relaxation, while smaller particles are considered to slip through the entanglement mesh. These theoretical results were confirmed by various experiments,20,23,29 while simulations based on bead−spring models were not able to capture such hopping.36 However, we have shown that direct accounting of the chemical structure of the polymer chains in the meltas implemented by fully atomistic modelsprovides a more detailed picture of nanoparticle diffusion. Specifically, in the present work we address diffusion in systems with ratio 2R/dT ∼ 0.1. For such systems, hopping is

still observed. Thus, in contrast to experiments concerning systems with 2R/dT ∼ 1.5−2, our work covers another limiting case, namely, the case of small 2R/dT ratios which could hardly be studied experimentally. Having established that hopping can be relevant both for systems with small 2R/dT ∼ 0.1 (the present work) and for large 2R/dT ∼ 1.5−2 (experiments), an important question now is how to extend the existing theoretical frameworks14,17 in order to strengthen the agreement between simulations, theory, and experiments. In view of our results, we believe that hopping should be considered as a general mechanism of diffusion in polymer melts for nanoparticles of sizes lying in the range between the typical correlation length and the reptation tube diameter. In summary, our simulations suggest that the diffusion of fullerenes is controlled by two mechanisms, i.e., segmental mobility of the chains at intermediate times and long-range hopping at longer times. Furthermore, our results suggest that hopping is relevant not only for particles whose diffusive properties are affected by the relaxation of the polymer entanglement mesh14,17 but also for nanoparticles whose dimensions are much smaller than the characteristic entanglement length. In other words, hopping diffusion may arise not only from the large-scale fluctuations of the mesh, but also from local interactions at the atomic level that can be taken into account by our all-atom simulations. In our view, this emphasizes the need to complement theoretical approaches with computational approaches able to provide a high level of detail of the relevant interactions. Finally, we showed that the FLE provides a reasonably good description of the ballistic and of the subdiffusive regime displayed by the motion of the fullerenes. While a firstprinciples approach for the normal diffusive regime is still lacking, a possible way to incorporate this in the framework of a phenomenological approach would be the inclusion of a superposition of exponential tails with different cutoff times in the kernel of the corresponding generalized Langevin equation. Developing a more comprehensive theoretical framework for the normal diffusion regime is beyond the scope of this paper, but it will be a major goal of our future research.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (S.V.L.). ORCID

Sergey V. Lyulin: 0000-0002-3743-4457 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Prof. S. B. Yuste for his careful reading of the manuscript as well as for helpful comments and suggestions. This study was supported by the Russian Foundation for Basic Research (RFBR), Grant No. 15-03-07614. The simulations were performed with the help of the computational resources of the Institute of Macromolecular Compounds of the Russian Academy of Sciences and the Lomonosov supercomputer at the Moscow State University. E.A. acknowledges financial support by Ministerio de Economiá y Competitividad (Spain) through Grant FIS2013-42840-P and by Junta de Extremadura through Grant GR15104. I


Article

Macromolecules

■

Nanoparticle Mobility in Semidilute Polymer Solutions. Macromolecules 2015, 48, 7585−7591. (25) Poling-Skutvik, R.; Krishnamoorti, R.; Conrad, J. C. SizeDependent Dynamics of Nanoparticles in Unentangled Polyelectrolyte Solutions. ACS Macro Lett. 2015, 4, 1169−1173. (26) Jee, A.-Y.; Curtis-Fisk, J. L.; Granick, S. Nanoparticle Diffusion in Methycellulose Thermoreversible Association Polymer. Macromolecules 2014, 47, 5793−5797. (27) Cherdhirankorn, T.; Best, A.; Koynov, K.; Peneva, K.; Muellen, K.; Fytas, G. Diffusion in polymer solutions studied by fluorescence correlation spectroscopy. J. Phys. Chem. B 2009, 113, 3355−3359. (28) Tuteja, A.; Mackay, M. E.; Narayanan, S.; Asokan, S.; Wong, M. S. Breakdown of the continuum stokes-einstein relation for nanoparticle diffusion. Nano Lett. 2007, 7, 1276−1281. (29) Grabowski, C. A.; Mukhopadhyay, A. Size Effect of Nanoparticle Diffusion in a Polymer Melt. Macromolecules 2014, 47, 7238−7242. (30) Grabowski, C.; Adhikary, B.; Mukhopadhyay, A. Dynamics of gold nanoparticles in a polymer melt. Appl. Phys. Lett. 2009, 94, 021903. (31) Mangal, R.; Srivastava, S.; Narayanan, S.; Archer, L. A. SizeDependent Particle Dynamics in Entangled Polymer Nanocomposites. Langmuir 2016, 32, 596−603. (32) Kalwarczyk, T.; Sozanski, K.; Ochab-Marcinek, A.; Szymanski, J.; Tabaka, M.; Hou, S.; Holyst, R. Motion of nanoprobes in complex liquids within the framework of the length-scale dependent viscosity model. Adv. Colloid Interface Sci. 2015, 223, 55−63. (33) Ganesan, V.; Pryamitsyn, V.; Surve, M.; Narayanan, B. Noncontinuum effects in nanoparticle dynamics in polymers. J. Chem. Phys. 2006, 124, 221102. (34) Liu, J.; Cao, D. P.; Zhang, L. Q. Molecular Dynamics Study on Nanoparticle Diffusion in Polymer Melts: A Test of the Stokes− Einstein Law. J. Phys. Chem. C 2008, 112, 6653−6661. (35) Pryamitsyn, V.; Ganesan, V. Dynamics of Probe Diffusion in Rod Solutions. Phys. Rev. Lett. 2008, 100, 128302. (36) Kalathi, J. T.; Yamamoto, U.; Schweizer, K. S.; Grest, G. S.; Kumar, S. K. Nanoparticle Diffusion in Polymer Nanocomposites. Phys. Rev. Lett. 2014, 112, 108301. (37) Starr, F. W.; Schrøder, T. B.; Glotzer, S. C. Effects of a nanoscopic filler on the structure and dynamics of a simulated polymer melt and the relationship to ultrathin films. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64, 021802. (38) Starr, F. W.; Schrøder, T. B.; Glotzer, S. C. Molecular Dynamics Simulation of a Polymer Melt with a Nanoscopic Particle. Macromolecules 2002, 35, 4481−4492. (39) Smith, J. S.; Bedrov, D.; Smith, G. D. A molecular dynamics simulation study of nanoparticle interactions in a model polymernanoparticle composite. Compos. Sci. Technol. 2003, 63, 1599−1605. (40) Kairn, T.; Daivis, P. J.; Ivanov, I.; Bhattacharya, S. N. Moleculardynamics simulation of model polymer nanocomposite rheology and comparison with experiment. J. Chem. Phys. 2005, 123, 194905. (41) Desai, T.; Keblinski, P.; Kumar, S. K. Molecular dynamics simulations of polymer transport in nanocomposites. J. Chem. Phys. 2005, 122, 134910. (42) Karatrantos, A.; Clarke, N.; Composto, R. J.; Winey, K. I. Polymer conformations in polymer nanocomposites containing spherical nanoparticles. Soft Matter 2015, 11, 382−388. (43) Li, Y.; Kröger, M.; Liu, W. K. Nanoparticle Effect on the Dynamics of Polymer Chains and Their Entanglement Network. Phys. Rev. Lett. 2012, 109, 118001. (44) Li, Y.; Kröger, M.; Liu, W. K. Dynamic structure of unentangled polymer chains in the vicinity of non-attractive nanoparticles. Soft Matter 2014, 10, 1723−1737. (45) Patti, A. Molecular Dynamics of Spherical Nanoparticles in Dense Polymer Melts. J. Phys. Chem. B 2014, 118, 3731−3742. (46) Patra, T. K.; Singh, J. K. Coarse-grain molecular dynamics simulations of nanoparticle-polymer melt: Dispersion vs. agglomeration. J. Chem. Phys. 2013, 138, 144901.

REFERENCES

(1) Bessonov, M. I.; Koton, M. M.; Kydryavtsev, V. V.; Laius, L. A. Polyimides: Thermally Stable Polymers; Plenum: New York, 1987. (2) Ghosh, M. Polyimides: Fundamentals and Applications; CRC Press: 1996. (3) Baker, A.; Dutton, S.; Kelly, D. Composite Materials for Aircraft Structures, 2nd ed.; AIAA: 2004. (4) Mouritz, A. P.; Gibson, A. G. Fire Properties of Polymer Composite Materials; Springer: Dordrecht, 2006. (5) Gay, D.; Hoa, S. V.; Tsai, S. W. Composite materials: Design and Application; CRC Press: London, 2003. (6) Rahmat, M.; Hubert, P. Carbon nanotube−polymer interactions in nanocomposites: A review. Compos. Sci. Technol. 2011, 72, 72−84. (7) Yoonessi, M.; Shi, Y.; Scheiman, D. A.; Lebron-Colon, M.; Tigelaar, D. M.; Weiss, R. A.; Meador, M. A. Graphene Polyimide Nanocomposites; Thermal, Mechanical, and High-Temperature Shape Memory Effects. ACS Nano 2012, 6, 7644−7655. (8) Kamanina, N. V.; Serov, S. V.; Shurpo, N. A.; Likhomanova, S. V.; Timonin, D. N.; Kuzhakov, P. V.; Rozhkova, N. N.; Kityk, I. V.; Plucinski, K. J.; Uskokovic, D. P. Polyimide-fullerene nanostructured materials for nonlinear optics and solar energy applications. J. Mater. Sci.: Mater. Electron. 2012, 23, 1538−1542. (9) Ajayan, P. M.; Tour, J. M. Materials Science: Nanotube composites. Nature 2007, 447, 1066−1068. (10) Wang, C.; Guo, Z.-X.; Fu, S.; Wu, W.; Zhu, D. Polymers containing fullerene or carbon nanotube structures. Prog. Polym. Sci. 2004, 29, 1079−1141. (11) Egorov, S. A. Anomalous nanoparticle diffusion in polymer solutions and melts: A mode-coupling theory study. J. Chem. Phys. 2011, 134, 084903. (12) Yamamoto, U.; Schweizer, K. S. Theory of nanoparticle diffusion in unentangled and entangled polymer melts. J. Chem. Phys. 2011, 135, 224902. (13) Yamamoto, U.; Schweizer, K. S. Microscopic Theory of the Long-Time Diffusivity and Intermediate-Time Anomalous Transport of a Nanoparticle in Polymer Melts. Macromolecules 2015, 48, 152− 163. (14) Dell, Z. E.; Schweizer, K. S. Theory of Localization and Activated Hopping of Nanoparticles in Cross-Linked Networks and Entangled Polymer Melts. Macromolecules 2014, 47, 405−414. (15) Brochard Wyart, F.; de Gennes, P. G. Viscosity at small scales in polymer melts. Eur. Phys. J. E: Soft Matter Biol. Phys. 2000, 1, 93−97. (16) Cai, L.-H.; Panyukov, S.; Rubinstein, M. Mobility of Nonsticky Nanoparticles in Polymer Liquids. Macromolecules 2011, 44, 7853− 7863. (17) Cai, L.-H.; Panyukov, S.; Rubinstein, M. Hopping Diffusion of Nanoparticles in Polymer Matrices. Macromolecules 2015, 48, 847− 862. (18) Omari, R. A.; Aneese, A. M.; Grabowski, C. A.; Mukhopadhyay, A. J. Diffusion of nanoparticles in semidilute and entangled polymer solutions. J. Phys. Chem. B 2009, 113, 8449−8452. (19) Guo, H.; Bourret, G.; Lennox, R. B.; Sutton, M.; Harden, J. L.; Leheny, R. L. Entanglement-Controlled Subdiffusion of Nanoparticles within Concentrated Polymer Solutions. Phys. Rev. Lett. 2012, 109, 055901. (20) Kohli, I.; Mukhopadhyay, A. Diffusion of Nanoparticles in Semidilute Polymer Solutions: Effect of Different Length Scales. Macromolecules 2012, 45, 6143−6149. (21) Kohli, I.; Alam, S.; Patel, B.; Mukhopadhyay, A. Interaction and diffusion of gold nanoparticles in bovine serum albumin solutions. Appl. Phys. Lett. 2013, 102, 203705. (22) Khorasani, F. B.; Poling-Skutvik, R.; Krishnamoorti, R.; Conrad, J. C. Mobility of Nanoparticles in Semidilute Polyelectrolyte Solutions. Macromolecules 2014, 47, 5328−5333. (23) Xue, C.; Zheng, X.; Chen, K.; Tian, Y.; Hu, G. Probing NonGaussianity in Confined Diffusion of Nanoparticles. J. Phys. Chem. Lett. 2016, 7, 514−519. (24) de Kort, D. W.; Rombouts, W. H.; Hoeben, F. J. M.; Janssen, H. M.; Van As, H.; van Duynhoven, J. P. M. Scaling Behavior of Dendritic J


Article

Macromolecules

translocation through lipid membranes. Nat. Nanotechnol. 2008, 3, 363−368. (67) Gaspar, H.; Fernandes, L.; Pereira, P.; Bernardo, G. Enhanced thermal stability of poly(methyl methacrylate) composites with fullerenes. Polym. Bull. 2015, 72, 1775−1786. (68) Pozdnyakov, A. O.; Kudryavtsev, V. V.; Friedrich, K. Sliding Wear of Polyimide-C60 Composite Coatings. Wear 2003, 254, 501− 513. (69) Kamanina, N. V.; Emandi, A.; Kajzar, F.; Attias, A.-J. LaserInduced Change in the Refractive Index in the Systems Based on Nanostructured Polyimide: Comparative Study with Other Photosensitive Structures. Mol. Cryst. Liq. Cryst. 2008, 486, 1043−1053. (70) Subocz, J.; Valozhyn, A.; Zenker, M. Effect of the carbon C60, C70 contents on the dielectric properties of polyimide films. Rev. Adv. Mater. Sci. 2007, 14, 193−196. (71) Yudin, V. E.; Svetlichnyi, V. M. Effect of the structure and shape of filler nanoparticles on the physical properties of polyimide composites. Russ. J. Gen. Chem. 2010, 80, 2157−2169. (72) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (73) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (74) Lyulin, S. V.; Larin, S. V.; Gurtovenko, A. A.; Lukasheva, N. V.; Yudin, V. E.; Svetlichnyi, V. M.; Lyulin, A. V. Effect of the SO2 group in the diamine fragment of polyimides on their structural, thermophysical, and mechanical properties. Polym. Sci., Ser. A 2012, 54, 631−643. (75) Lyulin, S. V.; Gurtovenko, A. A.; Larin, S. V.; Nazarychev, V. M.; Lyulin, A. V. Microsecond Atomic-Scale Molecular Dynamics Simulations of Polyimides. Macromolecules 2013, 46, 6357−6363. (76) Nazarychev, V. M.; Larin, S. V.; Lukasheva, N. V.; Glova, A. D.; Lyulin, S. V. Evaluation of the characteristic equilibration times of bulk polyimides via full-atomic computer simulation. Polym. Sci., Ser. A 2013, 55, 570−576. (77) Lyulin, S. V.; Larin, S. V.; Gurtovenko, A. A.; Nazarychev, V. M.; Falkovich, S. G.; Yudin, V. E.; Svetlichnyi, V. M.; Gofman, I. V.; Lyulin, A. V. Thermal properties of bulk polyimides: insights from computer modeling versus experiment. Soft Matter 2014, 10, 1224−1232. (78) Nazarychev, V. M.; Lyulin, A. V.; Larin, S. V.; Gurtovenko, A. A.; Kenny, J. M.; Lyulin, S. V. Molecular dynamics simulations of uniaxial deformation of thermoplastic polyimides. Soft Matter 2016, 12, 3972−3981. (79) Fal’kovich, S. G.; Larin, S. V.; Nazarychev, V. M.; Volgin, I. V.; Gurtovenko, A. A.; Lyulin, A. V.; Lyulin, S. V. Computer simulation of the heat-resistant polyimides ULTEM and EXTEM with the use of GROMOS53a6 and AMBER99 force fields. Polym. Sci., Ser. A 2014, 56 (4), 558−567. (80) Larin, S. V.; Falkovich, S. G.; Nazarychev, V. M.; Gurtovenko, A. A.; Lyulin, A. V.; Lyulin, S. V. Molecular-dynamics simulation of polyimide matrix pre-crystallization near the surface of a single-walled carbon nanotube. RSC Adv. 2014, 4, 830−844. (81) Larin, S. V.; Glova, A. D.; Serebryakov, E. B.; Nazarychev, V. M.; Kenny, J. M.; Lyulin, S. V. Influence of the carbon nanotube surface modification on the microstructure of thermoplastic binders. RSC Adv. 2015, 5, 51621−51630. (82) Falkovich, S. G.; Larin, S. V.; Lyulin, A. V.; Yudin, V. E.; Kenny, J. M.; Lyulin, S. V. Influence of the carbon nanofiller surface curvature on the initiation of crystallization in thermoplastic polymers. RSC Adv. 2014, 4, 48606−48612. (83) Falkovich, S. G.; Nazarychev, V. M.; Larin, S. V.; Kenny, J. M.; Lyulin, S. V. Mechanical Properties of a Polymer at the Interface Structurally Ordered by Graphene. J. Phys. Chem. C 2016, 120, 6771− 6777. (84) Flory, P. J. Statistical Mechanics of Chain Molecules; John Wiley and Sons: New York, 1969.

(47) Kalathi, J. T.; Grest, G. S.; Kumar, S. K. Universal Viscosity Behavior of Polymer Nanocomposites. Phys. Rev. Lett. 2012, 109, 198301. (48) Liu, J.; Gao, Y.; Cao, D.; Zhang, L.; Guo, Z. Nanoparticle Dispersion and Aggregation in Polymer Nanocomposites: Insights from Molecular Dynamics Simulation. Langmuir 2011, 27, 7926− 7933. (49) Cao, X.-Z.; Merlitz, H.; Wu, C.-X.; Ungar, G.; Sommer, J.-U. A theoretical study of dispersion-to-aggregation of nanoparticles in adsorbing polymers using molecular dynamics simulations. Nanoscale 2016, 8, 6964−6968. (50) Neyertz, S. Molecular Dynamics Simulations of Helium Permeation in Polyimides with a Bulky Dianhydride and a Fluorinated Diamine. Macromol. Theory Simul. 2007, 16, 513−524. (51) Hofman, D.; Ulbrich, J.; Fritsch, D.; Paul, D. Molecular modelling simulation of gas transport in amorphous polyimide and poly(amide imide) membrane materials. Polymer 1996, 37, 4773− 4785. (52) Chang, K.-S.; Tung, C.-C.; Wang, K.-S.; Tung, K.-L. Free Volume Analysis and Gas Transport Mechanisms of Aromatic Polyimide Membranes: A Molecular Simulation Study. J. Phys. Chem. B 2009, 113, 9821−9830. (53) Heuchel, M.; Hofmann, D.; Pullumbi, P. Molecular Modeling of Small-Molecule Permeation in Polyimides and Its Correlation to FreeVolume Distributions. Macromolecules 2004, 37, 201−214. (54) Kruse, J.; Kanzow, J.; Raltzke, K.; Faupel, F.; Heuchel, M.; Frahn, J.; Hofmann, D. Free Volume in Polyimides: Positron Annihilation Experiments and Molecular Modeling. Macromolecules 2005, 38, 9638−9643. (55) Neyertz, S.; Brown, D.; Pandiyan, S.; van der Vegt, N. F. A. Carbon Dioxide Diffusion and Plasticization in Fluorinated Polyimides. Macromolecules 2010, 43, 7813−7827. (56) Neyertz, S.; Brown, D. Influence of System Size in Molecular Dynamics Simulations of Gas Permeation in Glassy Polymers. Macromolecules 2004, 37, 10109−10122. (57) Ronova, I. A.; Rozhkov, E. M.; Alentiev, A. Y.; Yampolskii, Y. P. Occupied and Accessible Volumes in Glassy Polymers and Their Relationship with Gas Permeation Parameters. Macromol. Theory Simul. 2003, 12, 425−439. (58) Zhukova, E. K.; Kuznetsov, A. A.; Yablokova, M. Yu.; Alentiev, A. Yu. Gas separation properties of new thermoplastic polyimides with phenylamide groups in diamine moiety: Effect of polymer structure. Pet. Chem. 2014, 54 (7), 544−550. (59) Vidyakin, M. N.; Lazareva, Y. N.; Yampol’skii, Y. P.; Alentiev, A. Y.; Ronova, I. A.; Bruma, M.; Hamciuc, E.; Lungu, R. Transport properties of polyimides containing phenylquinoxaline moieties. Polym. Sci., Ser. A 2007, 49, 1045−1052. (60) Hanson, B.; Pryamitsyn, V.; Ganesan, V. Computer simulations of gas diffusion in polystyrene-C60 fullerene nanocomposites using trajectory extending kinetic Monte Carlo method. J. Phys. Chem. B 2012, 116, 95−103. (61) Vogiatzis, G. G.; Theodorou, D. N. Local Segmental Dynamics and Stresses in Polystyrene−C60 Mixtures. Macromolecules 2014, 47, 387−404. (62) Cao, Z.; Peng, Y.; Li, S.; Liu, L.; Yan, T. Molecular Dynamics Simulation of Fullerene C60 in Ethanol Solution. J. Phys. Chem. C 2009, 113, 3096−3104. (63) Bedrov, D.; Smith, G. D.; Davande, H.; Li, L. Passive transport of C60 fullerenes through a lipid membrane: a molecular dynamics simulation study. J. Phys. Chem. B 2008, 112, 2078−2084. (64) Li, L.; Davande, H.; Bedrov, D.; Smith, G. D. A Molecular Dynamics Simulation Study of C60 Fullerenes Inside a Dimyristoylphosphatidylcholine Lipid Bilayer. J. Phys. Chem. B 2007, 111, 4067− 4072. (65) Qiao, R.; Roberts, A. P.; Mount, A. S.; Klaine, S. J.; Ke, P. C. Translocation of C60 and its derivatives across a lipid bilayer. Nano Lett. 2007, 7, 614−619. (66) Wong-Ekkabut, J.; Baoukina, S.; Triampo, W.; Tang, I.-M.; Tieleman, D. P.; Monticelli, L. Computer simulation study of fullerene K


Article

Macromolecules (85) Hofling, F.; Franosch, T.; Frey, E. Localization Transition of the Three-Dimensional Lorentz Model and Continuum Percolation. Phys. Rev. Lett. 2006, 96, 165901. (86) Metzler, R.; Jeon, J.-H.; Cherstvy, A. G.; Barkai, E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys. 2014, 16, 24128−24164. (87) Sokolov, I. M. Models of anomalous diffusion in crowded environments. Soft Matter 2012, 8, 9043−9052. (88) Meroz, Y.; Sokolov, I. M. A toolbox for determining subdiffusive mechanisms. Phys. Rep. 2015, 573, 1−29. (89) Magdziarz, M.; Weron, A.; Burnecki, K.; Klafter, J. Fractional Brownian Motion Versus the Continuous-Time Random Walk: A Simple Test for Subdiffusive Dynamics. Phys. Rev. Lett. 2009, 103, 180602. (90) Metzler, R.; Jeon, J.-H.; Cherstvy, A. G. Non-Brownian diffusion in lipid membranes: Experiments and simulations. Biochim. Biophys. Acta, Biomembr. 2016, 1858, 2451−2467. (91) Tejedor, V.; Bénichou, O.; Voituriez, R.; Jungmann, R.; Simmel, F.; Selhuber-Unkel, C.; Oddershede, L. B.; Metzler, R. Quantitative Analysis of Single Particle Trajectories: Mean Maximal Excursion Method. Biophys. J. 2010, 98, 1364−1372. (92) Jeon, J.-H.; Metzler, R. Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries. Phys. Rev. E 2010, 81, 021103. (93) Jeon, J. H.; Leijnse, N.; Oddershede, L. B.; Metzler, R. Anomalous diffusion and power-law relaxation of the time averaged mean squared displacement in worm-like micellar solutions. New J. Phys. 2013, 15, 045011. (94) Sandev, T.; Tomovski, Z.; Dubbeldam, J. L. A. Generalized Langevin equation with a three parameter Mittag-Leffler noise. Phys. A 2011, 390, 3627−3636. (95) Goychuk, I. Viscoelastic subdiffusion: From anomalous to normal. Phys. Rev. E 2009, 80, 046125. (96) Kneller, G. R.; Baczynski, K.; Pasenkiewicz-Gierula, M. Communication: Consistent picture of lateral subdiffusion in lipid bilayers: Molecular dynamics simulation and exact results. J. Chem. Phys. 2011, 135, 141105. (97) Shin, H. K.; Kim, C.; Talkner, P.; Lee, E. K. Brownian motion from molecular dynamics. Chem. Phys. 2010, 375, 316−326. (98) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids, 4th ed.; Academic Press Inc.: New York, 2013. (99) Busselez, R.; Cerclier, C. V.; Ndao, M.; Ghoufi, A.; Lefort, R.; Morineau, D. . Discotic columnar liquid crystal studied in the bulk and nanoconfined states by molecular dynamics simulation. J. Chem. Phys. 2014, 141, 134902. (100) Boshoff, J. H. D.; Lobo, R. F.; Wagner, N. J. Influence of Polymer Motion, Topology and Simulation Size on Penetrant Diffusion in Amorphous, Glassy Polymers: Diffusion of Helium in Polypropylene. Macromolecules 2001, 34, 6107−6116. (101) Habasaki, J.; Ngai, K. L.; Hiwatari, Y. Molecular dynamics study of cage decay, near constant loss, and crossover to cooperative ion hopping in lithium metasilicate. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2002, 66, 021205. (102) Mukherjee, B.; Peter, C.; Kremer, K. Dual translocation pathways in smectic liquid crystals facilitated by molecular flexibility. Phys. Rev. E 2013, 88, 010502. (103) Takeuchi, H.; Okazaki, K. Dynamics of Small Molecules in a Dense Polymer Matrix: Molecular Dynamics Studies. Mol. Simul. 1996, 16, 59−74. (104) Kim, J.; Kim, C.; Sung, B. J. Simulation Study of Seemingly Fickian but Heterogeneous Dynamics of Two Dimensional Colloids. Phys. Rev. Lett. 2013, 110, 047801. (105) Patti, A.; El Masri, D.; van Roij, R.; Dijkstra, M. Stringlike Clusters and Cooperative Interlayer Permeation in Smectic Liquid Crystals Formed by Colloidal Rods. Phys. Rev. Lett. 2009, 103, 248304. (106) Lettinga, M. P.; Grelet, E. Self-Diffusion of Rodlike Viruses through Smectic Layers. Phys. Rev. Lett. 2007, 99, 197802.

(107) Wang, B.; Kuo, J.; Bae, S. C.; Granick, S. When Brownian diffusion is not Gaussian. Nat. Mater. 2012, 11, 481−485. (108) Chubynsky, M. V.; Slater, G. W. Diffusing Diffusivity: A Model for Anomalous, yet Brownian, Diffusion. Phys. Rev. Lett. 2014, 113, 098302. (109) Phillies, G. D. In complex fluids the Gaussian Diffusion Approximation is generally invalid. Soft Matter 2015, 11, 580−586. (110) Chaudhuri, P.; Berthier, L.; Kob, W. Universal Nature of Particle Displacements close to Glass and Jamming Transitions. Phys. Rev. Lett. 2007, 99, 060604. (111) Laage, D.; Hynes, J. T. On the Residence Time for Water in a Solute Hydration Shell: Application to Aqueous Halide Solutions. J. Phys. Chem. B 2008, 112, 7697−7701. (112) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987.

L


Molecular Dynamics Simulations of Fullerene Diffusion in Polymer

Recommend Documents