Binary Interactions between Dendrimer Molecules. A Simulation Study

Aug 1, 2014 - Departamento de Ciencias y Técnicas Fisicoquímicas, Facultad de Ciencias, Universidad Nacional de Educación a Distancia. (UNED), Pase...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/Macromolecules

Binary Interactions between Dendrimer Molecules. A Simulation Study Ana M. Rubio,† Carl McBride,‡ and Juan J. Freire*,‡ †

Departamento de Química Física, Facultad de Ciencias Químicas, Universidad Complutense, 28040 Madrid, Spain Departamento de Ciencias y Técnicas Fisicoquímicas, Facultad de Ciencias, Universidad Nacional de Educación a Distancia (UNED), Paseo Senda del Rey 9, 28040 Madrid, Spain



ABSTRACT: The effective interaction between two dendrimer molecules has been calculated from pairs of conformations, generated by Monte Carlo simulations, for coarse-grained models of PPI-DAB and PAMAM-EDA of different generation numbers (g = 2 to 7). The resulting interaction curves are compared with a theoretical Gaussian function. It is shown that the different simulation data can be adequately fitted to Gaussian functions in the relevant range of overlapping distances, but the fitting parameters are significantly different from those predicted by the theory for the highest values of g (which exhibit stronger interactions). Assuming a solvent mediated potential, the theoretical and Monte Carlo generated interaction curves are employed in molecular dynamics calculations of structure factors for different PPI-DAB dendrimer solutions having different concentrations. These calculations demonstrate that the different theoretical and simulation interactions may lead to significantly different quantitative and qualitative predictions for structural properties of dendrimer systems. Larger dendrimer molecules may be adequately described according to the observed trends in the parameter variations with generation number for these types of dendrimers.

I. INTRODUCTION The study of conformational properties of dendrimers has attracted a great deal of attention given due to potential application in various fields. Dendrimers are polymers showing a treelike branching with radial symmetry and multifunctional groups at their ends.1 The terminal groups can react to produce newer generations with a higher number of monomer units. Therefore, the generation number (g) which is defined by setting the value g = 1 for the core molecule (alternatively it is also described as G, with G = 0 for the core) is closely related to the molecule molecular weight, M. Most dendrimers are significantly smaller than regular polymers and exhibit a denser disposition of units. Numerical simulations provide a direct understanding of some of the dendrimer properties. In particular, they can lead to a deeper insight as to the similarities and differences for a variety of dendrimer types having diverse topologies, chemical compositions and, in the case of dendrimers with highly charged polar groups, pH or salt concentration.2,3 With this in mind, a great variety of molecular dynamics (MD) simulations with atomistic4,5 and united atom6 models, together with Brownian dynamics7,8 and Monte Carlo (MC) simulations with coarse-grained or bead models,9,10 have been undertaken. Most of the conformational studies have been directed to the study of an isolated dendrimer molecule. The effective interaction between dendrimer molecules has also been addressed from a theoretical point of view. This property is an important first step in understanding the properties of concentrated solutions. A theoretical description of the effective © 2014 American Chemical Society

interaction potential, weff (R), between two dendrimer molecules whose centers of masses are at a distance R has been proposed in terms of a Gaussian function11,12 ⎛ ⎞3/2 3 ⎟ exp( − 3R2/4R g . eff 2) weff (R )/kBT = Nt v0⎜⎜ 2⎟ 4 R π ⎝ g , eff ⎠ 2

(1)

where kBT is the Bolztmann factor, Nt the number of monomeric units, Rg,ef f 2 is the mean quadratic radius of gyration that corresponds to this description and v0 is a parameter describing the strength of the excluded volume for the monomeric units. This form assumes the Guinier theoretical behavior at low values of the scattering variable, leading to a convolution of two Gaussian profiles, and it is consistent with experimental data of the scattering structure factor for different concentrations.13,14 The shape of weff(R) for a pair of cluster-forming amphiphilic dendrimers placed at different distances between their centers of masses has also been explored via numerical simulation using a simple bead model.15 However, some characteristics of the dendrimer chemical structure (particularly, the degree of flexibility of the repeat units or linkers between branching points and the detailed description of the partial or total charges associated with the polar groups) may have a significant Received: May 30, 2014 Revised: July 9, 2014 Published: August 1, 2014 5379

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

procedure that we have previously employed in similar types of calculations for alkanes,21,22 two center Lennard-Jones models23 and polymers.24 We divide the MC trajectory in two sets of conformational samples and define the kth pair of conformations by selecting the kth conformation of each sample. For each pair of MC generated dendrimer conformations we locate the center of dendrimer 2 on a point on the surface of a sphere placed on the center of mass of dendrimer 1, with radius R, the fixed distance between the centers of masses. We choose the polar coordinates defining in terms of two random values for cos θ and Φ. Then, the intermolecular interactions between the two chains are obtained as

influence on the conformational behavior of single dendrimer molecules. These features have been effectively incorporated into the bead models. Thus, weff(R) has also been obtained by calculating the interacting forces between monomers of bead models of neutral or charged dendrimers.16 Nevertheless, it seems desirable to study the effective interaction with a model that can mimic the structure of particular real dendrimer molecules when introducing realistic parameters derived from atomistic models. In previous work we have investigated a coarse-grained model for MC simulations of dendrimers.10 This model consists of N hard-sphere beads located on the centers of the dendrimer branching units and in intermediate positions. The distribution of distances between these units has been introduced according to numerical data previously obtained from atomistic MD simulations for several types of dendrimer molecules. Following this method, we have obtained MC density profiles consistent with the MD simulations.17,18 Furthermore, we have obtained MC results for the radius of gyration, Rg, and intrinsic viscosities in good agreement with experimental data of polyamidoamine dendrimers with an ethylendiamine core (PAMAM-EDA) and polypropyleneimide with a diaminobutane core (PPI-DAB) in water as well as other types of dendrimers.10,19,20 In the present work, we calculate interactions for pairs of either PAMAM-EDA dendrimers or PPI-DAB dendrimers placed at different distances and orientations using conformations generated from our MC simulations. These two types of dendrimers and their derivatives are the most broadly studied and applied water-soluble dendrimers. We obtain conformational and orientational averages from which we can directly calculate the effective interaction potential and second virial coefficient for these dendrimers, without making any further assumptions as to the type of forces between the units. We calculate second virial coefficients which are compared with existing experimental data for the PPI-DAB case. Moreover, our data for the effective interaction potential are quantitatively compared with the theoretical Gaussian form. Both types of potentials have been subsequently used as solvent mediated point interactions in molecular dynamics simulations of PPIDAB dendrimer solutions at different concentrations. These simulations have allowed us to understand how the differences between theoretical and simulation interactions may affect experimentally accessible structural properties.

N

Uinter(R , θ , Φ) =

N

∑ ∑ U (R ij) i ,1

j ,2

(2)

where the summation terms correspond to values of the potential for intermolecular interactions between units i and j in dendrimers 1 and 2, respectively. U(Rij) adopts the hard-sphere values of infinity or zero for overlapping and nonoverlapping beads, Rij ≤ σ or Rij > σ. Parameter σ has been set at 4 and 5.8 Å for PPI-DAB and PAMAM-EDA10 respectively, in order to reproduce their experimental radius of gyration in water.25 We perform a certain number of evaluations of Uinter(R,θ,Φ) with different random angles, in order to fully explore the different orientations for every particular pair of conformations. From the results corresponding to the different pairs of conformations and orientations we evaluate the averaged Mayer function for the given value of R fM (R ) = ⟨exp( − Uinter(R , θ , Φ)/kBT )⟩conf, θ , Φ − 1

(3)

obtaining arithmetic means over all the sampled pairs of differently oriented conformations. We consider different values of R, obtaining the corresponding values of Uinter(R,θ,Φ) and perform further averages to calculate f M(R). Finally, we obtain the second virial coefficient as B2 = −2π

∫0



fM (R )R2 dR

(4)

Numerical integrations are truncated at a upper limit, R = Rmax beyond which the calculated Mayer function is consistently zero. It varies within the range Rmax = 60−160 Å for the g = 2− 7 dendrimers. We have performed 10 independent evaluations, or batches, of B2. Each batch uses 5−10 random orientational angles starting with a different choice of their initial values. Also, for each value of R, we obtain the effective interaction potential as

II. MODELS AND METHODS Parameters for the coarse-grained model corresponding to the PAMAM-EDA and PPI-DAB dendrimers have been provided in our previous work.10,19 The molecules are represented by hard-sphere beads placed on the branching atoms (nitrogen atoms of the amine groups) and at an intermediate point along each spacer. The hard-sphere bead diameter is chosen to provide a good reproduction of the dendrimer size for generations in the range g = 2−7 (or G = 1−6). The model uses a distribution of distances between beads obtained from previous atomistic MD simulations performed for smaller values of g. The trajectories extend over a total of 2 × 108 Monte Carlo steps, or move attempts (each move consists of a small displacement of single beads). We have saved a fraction (1/10 000) of these MC moves, i.e. a total of 2 × 104 conformations. We calculate the effective binary interaction potential and the second virial coefficient of dendrimers according to the

weff (R )/kBT = −ln[1 + fM (R )]

(5)

Furthermore, we use the theoretical or Monte Carlo generated effective interaction curves as solvent-mediated point potentials to simulate dendrimer solutions at different concentrations. Using the GROMACS software package,26 we have performed molecular dynamics simulations for systems composed of 1000 dendrimer particles. These molecules interact according to the assigned potentials in the NVT ensemble with cubic periodic boundary condition in boxes whose size are set to reproduce the desired concentrations. These simulations employ the Berendsen thermostat27 with lengths of 20−50 ns with a time step of 1 fs. Structure factors are obtained from the radial distribution functions calculated 5380

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

Results for the second virial coefficient can also be discussed in terms of the interpenetration factor,

from the dendrimer coordinates saved along the molecular dynamics trajectories, according to standard numerical procedures included in GROMACS.

Ψ* = 2B2 /(4πR g 2)3/2

III. RESULTS AND DISCUSSION In Figure 1, we show our simulation results obtained for the second virial coefficient, (in Å3) for PPI-DAB (σ = 4 Å),

(7)

that combines this magnitude with Rg2, the mean quadratic radius of gyration of the molecule.29 Table 1 shows the values obtained for Ψ* from the simulation results for B2 and the previously reported simulation Table 1. Values of Ψ* Calculated According to Equation 7 dendrimer PPI-DAB

PAMAM-EDA

Ψ*(simulation)

2 3 4 5 6 7 2 3 4 5 6 7

1.24 1.37 1.51 1.64 1.70 1.65 1.51 1.57 1.56 1.71 1.74 1.66

Ψ*(experimental) 1.2 1.2 1.0 2.0

± ± ± ±

0.2 0.2 0.2 0.4

Ψ* is obtained by combining simulation data(σ = 4 Å) in Figure 1 and experimental results28 (LS data for PPI-DAB) of the second virial coefficient and our simulation values for the radius of gyration.10 The relative errors of the theoretical results are always smaller than 5%. See text for further details.

Figure 1. Comparison between our simulation results(lines) and experimental data28(symbols) for the second virial coefficients of PPIDAB dendrimers (in Å3). Dashed, dot-dashed, solid, dash-dotted, and dotted lines correspond to σ = 2, 3, 4, 5, and 6 Å. Blue solid squares: LS experimental data, 25 °C. Red solid circles: VPO experimental data, 45 °C. Magenta open triangles: VPO data including a correction(see text). Simulation error bars are always smaller than 1% of the values.

results for Rg.10 They are compared with interpenetration factors obtained from the experimental results for the second virial coefficients28 (LS data, which exhibit smoother variations with g) combined with our simulation results for the radius of gyration for PPI-DAB. It can be observed that there are minor differences in the simulation values of Ψ* with g. However, a systematic variation, including a maximum, can be observed for both types of dendrimers. The experimental values show greater nonsystematic variations and higher uncertainties. It is interesting to note that some of the simulation and experimental results for Ψ* for the highest generation numbers are greater than the value for the hard-sphere model Ψ* = 1.62 (or, expressed in terms of the ratio between the second virial coefficient and the molecular volume of a compact spherical structure, B2/V = 4). For this theoretical model, molecules cannot penetrate at a distance from the sphere center smaller than its radius, rsph = (5/3)1/2Rg. In the case of a similarly impenetrable hollow sphere, rsph = Rg. This implies a radius of gyration greater than that of a hard sphere, leading to the smaller value of the parameter, Ψ* = 0.75. Linear flexible chains can easily interpenetrate each other because of their low total density of polymer units. Consequently, their second virial coefficients are notably smaller. The value Ψ* ≅ 0.26 was obtained for long flexible linear polymer coils in a good solvent.24 For flexible star polymers, parameter Ψ* increases with the number of arms (i.e., with increasing degree of compactness), approaching the hard sphere limit.24 However, the hard-sphere value should not be considered as an upper bound for Ψ*. For example, B2/V > 4 for hard ellipsoids.30 However, the asphericity (or deviation form the spherical shape)31 is small or moderate in the case of dendrimers, decreasing for higher generation numbers.32

compared with experimental data for different generations.28 These experimental results were obtained via two different techniques, low angle laser light scattering (LS) and vapor pressure osmometry (VPO), with or without a technical correction due to the presence of the methanolate ions, in methanol solutions (the experimental molecular weights for g = 4 and 5 samples are slightly smaller than the theoretical values, which may indicate uncomplete polymerization reactions for these cases.) The experimental data for the radius of gyration of these dendrimers28 are similar to those obtained from our Monte Carlo simulations,10 also consistent with PPI-DAB data measured in aqueous solutions.25 The experimental virial coeffcients, initially expressed as A2 (in units of cm3 mol/g2), are converted to B2 (in molecular volume units) via

B2 = A 2 M2 /NA

G

(6)

where NA is the Avogadro number and M is derived from the dendrimer chemical composition. It can be observed that the simulation results are in fair agreement with experimental data, with high uncertainties and nonsystematic variations between generations greater than g = 3. We have also plotted four more alternative simulation curves that correspond to the values σ = 2, 3, 5, and 6 Å, i.e. considering that σ can differ from the value employed for the intramolecular potential in the previous generation of the single dendrimer trajectories, in order to provide a more precise description of intermolecular interactions. The different curves are not far from the experimental data. Taking all the experimental data into consideration, σ = 4 Å seems a reasonable choice. 5381

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

Considering the cases described above, Ψ* is essentially determined by the behavior of the dendrimer density profile referred to the distance from the molecule center of masses, ρ(r). A greater density profile ρ(r) in the external regions of the dendrimer tends to increase the radius of gyration leading to a smaller Ψ* value, as in the case of impenetrable hollow spheres where ρ(r) is zero except at r = Rg, i.e., on the sphere surface. This is also the case of high molecular weight linear and star polymers in good solvents, which show decreasing density profiles,33 ρ(r) ∼ r‑4/3. Many dendrimer molecules show density profiles that may be tentatively described by a “dense core”, i.e. decreasing ρ(r), model.12,34 However, depending on their composition, the generation number, and other factors such as the pH or salt concentration in water5,18 they can exhibit hollow internal regions or a denser external shell. Therefore, ρ(r) may not be monotonous along the whole interval of r values, and in fact, it can be noticeably different from zero for relatively high values of r/Rg. This means that many dendrimer molecules are difficult to penetrate even at distances r/Rg greater than the (5/3)1/2Rg value of an equivalent hard-sphere. In Figure 2, we present the

Figure 3. ln[weff(R)/kBT] vs (R/Rg)2 for PPI-DAB dendrimers: (a) g = 2; (b) g = 7. Symbols: simulation results. Line: linear fit of initial slope.

always downward and they can be explained by taking into account that two dendrimers cannot interact when they are placed at distances far away from each other, while the Gaussian function shows a slow continuous decay. However, it is not clear if further interactions could occur between real atoms beyond the extension of the most external beads of our coursegrained model and they could contribute to extend the range for which the Gaussian behavior is obeyed for real molecules. The lower limit in the R/Rg range of Gaussian behavior marks the point where overlap between monomers occurs in all of the explored pairs of conformations. This limit increases with g, varying from 1.7 to 2.25 for the PAMAM-EDA dendrimers, a similar increase is also seen, from 1.1 to 2.15 for the PPI-DAB dendrimers. Therefore, we conclude that the range Gaussian behavior can vary significantly for different types of dendrimers or generation numbers. Thus, for the smallest PPI-DAB dendrimers, with less flexible and shorter spacers, the two dendrimers may show an important degree of overlap. Indeed, their centers of masses can be at a distance similar to the radius of gyration. However, dendrimers cannot approach a distance of more than twice Rg when the generation number increases. We have compared the performance of the theoretical Gaussian description of the effective interaction potential and other simple empirical functions that can be tentatively

Figure 2. Density profile, for the PAMAM-EDA g = 5 dendrimer. The mark corresponds to r = (5/3)1/2Rg.

density profile obtained from previous simulation data17 for our coarse-grained model of PAMAM-EDA with g = 5. It is shown that the density is significant at (5/3)1/2Rg and, therefore, overlap between the dendrimers still can be found at distances beyond this value. Therefore, a relatively high density at intermediate r values, which extends beyond the equivalent hard-sphere contact distance, can explain the higher values of Ψ* for dendrimers at moderate or high generation numbers. Using eq 5, we obtain the effective interaction potential curves from the simulation data corresponding to our different dendrimer molecules. lnweff(R) vs (R/Rg)2 representations are employed to check the theoretical Gaussian form provided by eq 1. Figures 3 and 4 show these representations for PPI-DAB and PAMAM-EDA molecules (low and high g). It can be observed that the linear fit is consistent with eq 1 along a welldefined range of R/Rg values. The upper limit of this range is placed at R/Rg =2.6−2.75, with its precise value following some slight variations depending on g and the dendrimer type. The smaller values correspond to the highest generation number. Also, PAMAM-EDA dendrimers show slightly earlier deviations from the Gaussian behavior than PPI-DAB. Deviations are 5382

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

possible empirical forms of the effective interaction potential is clearly demonstrated for the low g dendrimers. Meanwhile,, the high g dendrimers behave more like hard spheres for both types of molecules. As previously discussed,11,13,14 the dendrimer density profile changes drastically as the generation number grows and a Gaussian effective interaction cannot be obtained as the result of the theoretical convolution if the individual profiles are not Gaussian-like themselves. The numerical parameters obtained in our linear fits of the lnweff (R) vs (R/Rg)2 representations within the adequate ranges can be discussed according to the theoretical description provided by eq 1. Thus, the fitted slopes directly yield an estimate for the ratio Rg,eff/Rg. In Figure 5, we show the ratios

Figure 5. Rg,eff/Rg vs g for PPI-DAB(red circles) and PAMAMEDA(blue squares) dendrimers.

between these numerical estimations of this ratio between the radius of gyration from the effective potential and the values of Rg directly calculated from our simulation trajectories. The ratio is slightly greater than 1 for the lowest generations and decreases, becoming significantly smaller than 1 for the high g dendrimers. The highest differences between Rg,eff and Rg are about 25%. It should be noticed that these differences are directly related with deviations from the theoretical form of the effective interactions, described by eq 1. It can be observed that the ratios corresponding to both types of dendrimer show fairly similar trends and they are quantitatively close for most generation numbers, even though the chemical structures of PAMAM-EDA and PPI-DAB are different (linkers are substantially more flexible in the PAMAM-EDA case). This suggests the use of similar values of Rg,eff/Rg at different generations to describe the effective potential of these and other similar dendrimer structures. From ordinates of our linear fits for ln weff (R) vs (R/Rg)2, together with the numerical values of Rg,ef f, we have calculated the parameter representing the excluded volume of the monomeric units, ν0, according to eq 1. In the case of our model, these units are hard-sphere beads. Therefore, the bead diameters can be estimated from

Figure 4. ln[weff(R)/kBT] vs (R/Rg)2 for PAMAM-EDA dendrimers: (a) g = 2 and (b) g = 7. Symbols: simulation results. Line: linear fit of initial slope.

employed for the same purpose. Representations of the type lnweff(R) vs (R/Rg), associated with an exponential decay deviate sooner from the simulation points. The performance is not so good with respect to the Gaussian form, covering a similar range of values of R/Rg for the high g dendrimers, while the upper limit of an exponential function is significantly lowered, i.e. the interval is substantially narrowed, in the case of the smallest dendrimers (e.g., the upper limit of the exponential decay is R/Rg = 2 for both types of dendrimers with g = 2). Representations of lnweff(R) vs (R/Rg)n with n > 2 cannot be precisely fitted to a linear function. With n = 3, this problem can only be detected via the numerical analysis of the fitting correlations for the highest values of g but, in the case of the smallest dendrimers, moderate or important nonlinear curvatures can be visually appreciated for PPI-DAB and PAMAM-EDA. Consequently, we conclude that our simulation model shows an adequate agreement with the Gaussian behavior in the range of values of R corresponding to a slight or moderate repulsion between dendrimers. The description is better for the dendrimers with low g that can interpenetrate easily, especially in the case of the PPI-DAB. Moreover, better performance of the theoretical Gaussian description with respect to other

ν0 = (2π /3)σeff 3

(8)

and we can compare these estimates with the bead diameter values initially assigned in the simulations. In Figure 6 we show the ratios σeff/σ between bead diameters obtained from the 5383

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

Figure 7. Internal averaged density d (in macrosocopic units, g/cm3) vs g for PPI-DAB(red circles) and PAMAM-EDA(blue squares) dendrimers.

Figure 6. σeff/σ vs g for PPI-DAB(red circles) and PAMAM-EDA(blue squares) dendrimers.

force derived directly from the model parameters. The results are presented in Figures 8 and 9, where we show the weff(R)/

Gaussian interaction potential and those considered in the simulations for the PPI-DAB and PAMAM-EDA dendrimers. It can be observed that the greatest difference between both types of diameters is about 20% for the g = 7 PAMAM-EDA dendrimer. Therefore, there is a fair agreement between the model bead diameter and the value needed for a correct description of Gaussian interactions between units in all cases. However, there is a systematic variation of the ratios σeff/σ with both the generation number and the type of dendrimer. The qualitative aspect of the variation with g is similar for both types of dendrimers. Thus, the ratios are smaller than 1 for small g, have a minimum for the intermediate size dendrimers and are noticeable greater than 1 for the highest g values. The PPI-DAB ratios are always smaller than those corresponding to the PAMAM-EDA dendrimers. The change in the effective excluded volume may be related with the bead density, since a higher internal bead density implies a smaller number of contacts with external units. In fact, the qualitative trend observed in Figure 6 is similar to the variation of the averaged dendrimer internal bead density, d, obtained by assuming a compact sphere containing the mass of a dendrimer molecule, i.e. d= M/[4π(5/3)3/2NARg3/3]. These data are shown in Figure 7. However, it should be also noted that the PPI-DAB dendrimers with lower σeff/σ ratios have a greater density than the PAMAM-EDA dendrimers for a given generation number (even though their monomeric mass, not considered in the simulation model, is smaller). Consequently, a density correlation cannot explain the different ratios found for two different dendrimers with a given generation number. The observed correlation between effective bead diameters and densities cannot be clearly described from the quantitative point of view and it is only valid for each particular type of dendrimer. However, the similar qualitative variation with the generation number permits a common description of different dendrimers, using ν0 as an empirical parameter that can be fitted to describe experimental data for different dendrimersolvent systems.14 We will provide a practical example of this procedure below. The simulated and fitted values of weff(R) can also be more directly compared with the predictions obtained with eq 1 and eq 8 using the parameter values Rg,eff = Rg and σeff = σ values, wG(R). This way, we can check the accuracy of an effective

Figure 8. Ratios between simulation(symbols) or fitted Gaussian(lines with the same color codes) results for the effective interaction potential of the PPI-DAB dendrimers. The dashed line corresponds to a linear fit for the g = 2 case. Symbols: g = 2, dark red upward triangles; g = 3, red downward triangles; g = 4, orange diamonds; g = 5, green stars; g = 6, blue circles; g = 7, purple squares. A line marking the overlapping distance of a pair of equivalent spheres is also shown.

wG(R) ratios obtained for the different cases. These ratios show significant deviations from 1, with similar trends for both types of dendrimers. The deviations are greater for the highest generations that show weff(R) values substantially smaller than those predicted by the theory. In the same figures it can be verified that similar ratios are obtained when the simulation values are substituted by the fitted Gaussian function, i.e., calculated with our previously described values of Rg,eff and σeff, in the relevant overlapping ranges where they were fitted. Fitting ranges can be easily discerned and compared in these particular figures. A maximum is shown for the fitting functions of low generation dendrimers, marking the upper limit of the Gaussian range. Moreover, we can verify the significant reduction of the upper limit of fitting interval when a linear 5384

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

preformed with a different hard-sphere bead diameter, σ = 2.25 Å, which is consistent with our (ν0)corr value. Our results for the structure factors are obtained from molecular dynamics trajectories of point particles whose interactions are described by the Gaussian fits to our simulation results. (The technical characteristics of the molecular dynamics simulations are given above). The concentration for the different dendrimers are expressed as ρ, in beads (of our simulation model) per Å3. The chosen concentrations are referenced to the overlapping concentration for rigid spheres of the same radius of gyration, ρ*=N/[(4π/3)(5/3)3/2Rg3]. In Figure 10, we show the results corresponding to g = 4 for two different concentrations. The general aspect of our curves

Figure 9. As Figure 8 but for PAMAM-EDA dendrimers.

function is used instead of a Gaussian for the case of the g = 2 dendrimers discussed above (though the linear representation slightly improves the fit in the lower limit for the particular case of the PPI-DAB dendrimer). Taking into account all the different cases, we confirm that a Gaussian effective interaction potential is consistent with our simulation model, but only when adequately fitted parameters are used. In the case of the intermediate generation numbers, Rg,eff ≅ Rg implying that the difference between theoretical and simulation results only affects the pre-exponential factor. Structure factors obtained via Monte Carlo numerical simulations for systems of coarse-grained dendrimers, g = 4, at different concentrations have been previously performed, showing a good agreement with simulations obtained with the theoretical effective interactions between the centers of masses.35 Consequently, many-body effects seems to be small, indicating that the effective interactions may be able to give a good structural description of dendrimer systems. Moreover, good agreement has been previously demonstrated14 between structure factors, S(q), derived from the theoretical eq 1 for interaction between dendrimers (following a theoretical calculation based in the Percus−Yevick closure approximation) and experimental scattering data for different concentrations of particular functionalized g = 4 PPI-DAB dendrimer (with end groups −NH−C6H5) in the solvent dimethylacetamide (DMA). This comparison was carried out using appropriately chosen values of the solvent-mediated excluded volume parameter, (ν0)L, for the theoretical beads (93 Å3 for the highest concentration and 76 Å3 for the remaining cases) and using a theoretical model of Nt = 62 beads. Mimicking this approach, we have obtained structure for our PPI-DAB dendrimer model using our fitted interaction potentials but introducing a value of the excluded volume parameter (ν0)corr = 24 Å3 equivalent to (ν0)L = 93 Å3 for our N = 122 bead simulation model, N2(v0)corr = Nt2(v0)L. This leads to weaker interactions between dendrimers that can also be described by introducing a correction factor of fcorr = 0.336 in the fitted preexponential parameters in order to reproduce the theoretical value. An equivalent method has been previously employed to predict interaction curves by combining the profiles obtained from single dendrimer MC simulations with an empirical excluded volume parameter.36 As a second option, we have considered that the intermolecular interaction between these particular molecules can be described through simulations

Figure 10. Structure factor for the g = 4 dendrimer obtained from molecular dynamics simulations from fitted Gaussian function interactions, with the pre-exponential factors corrected with fcorr = 0.336(up triangles), or performed with σ = 2.25 Å (down triangles), and from the theoretical Gaussian interaction(lines) calculated with(ν0)corr = 24 Å3. Solid line and symbols (in red): ρ = 2.9 × 10−3 beads/Å3. Dashed line and open symbols(in blue): ρ = 8.7 × 10−4 beads/Å3. ρ* = 9.1 × 10−3 beads/Å3.

in Figure 10 do not differ much from the theoretical calculation of structure factors for the functionalized g = 4 PPI-DAB dendrimer reported previously14 using two values of ρ that correspond to the experimental measurements mentioned above for dendrimer volume fractions Φ = 0.237 and Φ = 0.070. (The standard g = 4 PPI-DAB dendrimer has a smaller size than the functionalized dendrimer and, therefore, both values of ρ are in our case are well below overlapping concentration). The agreement between the results obtained through eq 1 and from the interaction curve fitted from our simulation data with the pre-exponential correction factor, fcorr = 0.336, is very good. This should be expected since, in this particular case, the interaction function pre-exponential factors are identical by definition and Rg,eff ≅ Rg. It is verified that the curves obtained from the simulations performed with σ = 2.25 Å are similarly close. For other values of g, however, Rg,eff and Rg show significant differences. The pre-exponential factors can also be different since the ratios between theoretical and fitted function ratios vary from those obtained for g = 4. Consequently, we have investigated whether the theoretical interaction potentials calculated with the choices Rg,ef f = Rg and ν0 = (ν0)corr in eq 1 can also reproduce the structure factors obtained with the fitted functions (with pre-exponential terms modified by the factor 5385

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

fcorr = 0.336 or obtained directly from the simulations with σ = 2.25 Å) for the g = 7 and g = 2 PPI-DAB molecules. In Figure 11, we show the structure functions obtained for g = 7. We also consider two concentrations below the

Figure 12. Structure factor for the g = 2 dendrimer obtained from molecular dynamics simulations from fitted Gaussian function interactions, with the pre-exponential factors corrected with fcorr = 0.336 (up triangles), or performed with σ = 2.25 Å (down triangles), and from the theoretical Gaussian interaction (lines) calculated with (ν0)corr = 24 Å3. Solid line and symbols (in red): ρ = 2.6 × 10−2 beads/ Å3. Dashed line and open symbols (in blue): ρ = 3.2 × 10−3 beads/Å3. ρ* = 1.0 × 10−2 beads/Å3.

Figure 11. Structure factor for the g = 7 dendrimer obtained from molecular dynamics simulations from fitted Gaussian function interactions, with the pre-exponential factors corrected with fcorr = 0.336 (up triangles), or performed with σ = 2.25 Å (down triangles), and from the theoretical Gaussian interaction (lines) calculated with (ν0)corr = 24 Å3. Solid line and symbols (in red): ρ = 8.1 × 10−3 beads/ Å3. Dashed line and open symbols (in blue): ρ = 3.0 × 10−3 beads/Å3. ρ* = 1.2 × 10−2 beads/Å3.

remarkably different descriptions of the solvent mediated excluded volume effects. This suggests than a single simulation of binary interactions with a reasonable choice for the hardsphere diameter for intermolecular interactions is sufficient to describe the potential curves for solutions in different solvent through a Gaussian fit of the simulation data and adequate corrections of the pre-exponential parameter. It should be stressed that, as described above, the value (ν0)corr = 24 Å3, or (ν0)corr = 93 Å3 for the solvent mediated excluded volume is consistent with scattering experimental data for a particularly functionalized PPI-DAB dendrimer in DMA at a given concentration.14 The phenyl rings in the −NH−C6H5 end groups are mutually attractive favoring the approach between dendrimer molecules. Therefore, this excluded volume parameter may not be correct for the standard PPI-DAB dendrimer (with positively charged amine groups at the ends) in other polar solvents. Comparing the different simulation data in Figure 1, we can observe that the bead diameter σ = 2.25 Å curve, consistent with (ν0)corr = 93 Å3, would be located in the lowest bound of the range that could be considered compatible with the existing experimental data virial coefficients in methanol.28 Consequently, higher volumes of the excluded volume parameter and the bead diameter are probably more realistic. This would imply stronger intermolecular interactions and shaper structure factors. In summary, our simulations with a hard-sphere bead model that incorporate the short-range characteristics of real dendrimers and reproduce their experimental sizes10 show binary intermolecular interactions compatible with existing experimental second virial coefficient data of PPI-DAB molecules of different generation numbers.28 The resulting effective interaction potential for different PPI-DAB and PAMAM-EDA dendrimers can be interpreted in terms of a theoretically derived Gaussian form11 in the specific range of significant interacting distances between the dendrimer centers of masses for all generations. The extent of this range has a noticeable dependence on the generation number and type of

overlapping value. In both cases, the agreement between the results based on simulation data and the theoretical values is not so good. For the more concentrated system, dendrimers clearly behave as hard spheres with two prominent peaks in the considered range of q values. eq 1 is able to predict the location of the first peak though it overestimates the height and width obtained from the simulation results. The smaller effective radius of gyration derived form the simulations yields a faster decay of interactions, leading to weaker fitted interactions at overlapping distances from the simulation data. The more dilute system only shows a prominent peak compatible with a softer potential, but eq 1 overestimates its height and location, providing a poor description of the interactions. The limitations of the theoretical description for higher generation number dendrimers discussed above can explain these discrepancies. In the case of the g = 2 dendrimer, shown in Figure 12, the pre-exponential factor for eq 1 and from our fits to the simulation data (corrected with fcorr = 0.336 or obtained with σ = 2.25 Å) is small. In this case, the values of Rg,eff (derived from the simulation curves) and Rg are closer. Moreover, interactions are softer for these smaller dendrimers which can easily interpenetrate. Consequently, eq 1 gives structure factors that are remarkably close to those obtained with the two types of simulation data. The agreement is good even for the higher concentration, which was chosen to be well above the overlapping limit. It is interesting to note that the two alternative simulation curves (obtained with the pre-exponential correction factor fcorr = 0.336 or directly fitted from the simulation data corresponding to σ = 2.25 Å) are remarkably similar. We only note a marginal difference; the σ = 2.25 Å curves are flatter for the most diluted solvents implying a slightly higher osmotic compressibility. This agreement is found in spite of their 5386

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387

Macromolecules

Article

(17) Freire, J. J. Polymer. Soft Matter 2009, 5, 1912. (18) Freire, J. J.; Ahmadi, A.; McBride, C. J. Phys. Chem. B 2013, 117, 15157. (19) Rodriguez, E.; Freire, J. J.; del Rio Echenique, G.; Hernández Cifre, J. G.; García de la Torre, J. Polymer 2007, 48, 1155. (20) del Rio Echenique, G.; Rodríguez Schmidt, R.; Freire, J. J.; Hernández Cifre, J. G.; García de la Torre. J. Am. Chem. Soc. 2009, 131, 8548. (21) López Rodríguez, A.; Freire, J. J. Mol. Phys. 1988, 63, 591. (22) López Rodríguez, A.; Vega, C.; Freire, J. J.; Lago, S. Mol. Phys. 1991, 73, 691. (23) Vega, C.; McBride, C.; Menduiña, C. Phys. Chem. Chem. Phys. 2002, 4, 3000. (24) Rubio, A. M.; Freire, J. J. Macromolecules 1996, 29, 6946. (25) Scherrenberg, R.; Coussens, B.; van Vliet, P.; Edouard, G.; Brackman, J.; deBrabander, E.; Mortensen, K. Macromolecules 1998, 31, 456. (26) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (27) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon: Oxford, England, 1987. (28) Rietveld, I. B.; Smit, J. A. M. Macromolecules 1999, 32, 4608. (29) Yamakawa, H. Macromolecules 1992, 25, 1912. (30) McBride, C.; Lomba, E. Fluid Phase Equilib. 2007, 255, 37. (31) Wei, G.; Eichinger, B. E. J. Chem. Phys. 1990, 93, 1430. (32) Freire, J. J.; Rubio, A. M. Polymer 2008, 49, 2762. (33) Freire, J. J. Adv. Polym. Sci. 1999, 143, 35. (34) Boris, D.; Rubinstein, M. Macromolecules 1996, 29, 7251. (35) Götze, I. O.; Likos, C. N. J. Phys.: Condens. Matter 2005, 17, S1777. (36) Harreis, H. M.; Likos, C. N.; Ballauff, M. J. Chem. Phys. 2003, 118, 1979.

dendrimer, becoming narrower as the dendrimer size increases. These limitations can be understood as a consequence of the failure of the theoretical description in terms of Gaussian-like profiles. Furthermore, the effective values of the dendrimer radius of gyration derived by fitting to a Gaussian function show some deviations with respect to the results evaluated from our initial Monte Carlo simulations. These differences are similar for the two types of dendrimers when a given generation number is considered. Further differences are found between the effective bead diameter derived from the excluded volume parameter of the Gaussian function and those initially assigned in the model. The diameter is related to the excluded volume strength, which can be fitted to experimental data. Computing the structure factors derived from interactions functions obtained from theory or from the MC simulation data for particular PPI-DAB solutions, we have shown that there are significant differences for the dendrimers with high generation number. Therefore, adequate parameters should be taken in consideration in order to model an effective Gaussian intermolecular potential that can properly describe the interactions of the dendrimers. The present simulations show that the parameters, effective radius of gyration and bead diameter (or excluded volume strength), display sufficient similarity for these two broadly studied and applied dendrimers as to suggest that a general procedure can be implemented to describe them and other similar dendrimer systems.



AUTHOR INFORMATION

Corresponding Author

*E-mail: (J.J.F.) [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work has been partially supported by Research Project CTQ2010-16414 from the Spanish Government, DGICT, and a UNED Postdoctoral Grant, 2013-018-UNED-POST, for C.M.



REFERENCES

(1) Fréchet, J. M. J., Tomalia, D. A., Eds. Dendrimers and other dendritic polymers; Wiley: Chichester, England, 2001. (2) Likos, C. N.; Ballauff, M. Angew. Chem., Int. Ed. 2004, 43, 2998. (3) Freire, J. J. Soft Matter 2008, 4, 2139. (4) Karatasos, K. Macromolecules 2005, 38, 4472. (5) Maiti, P. K.; Ç ağin, T.; Lin, S.-T.; Goddard, W. A., III Macromolecules 2005, 38, 979. (6) Lyulin, S. V.; Evers, L. J.; van der Schoot, P.; Darinskii, A. A.; Lyulin, A. V.; Michels, M..A. J. Macromolecules 2004, 37, 3049. (7) Bosko, J. T.; Prakash, J. R. Macromolecules 2011, 44, 660. (8) Timoshenko, E. G.; Kuznetsov, Y. A.; Connolly, R. J. Chem. Phys. 2002, 117, 9050. (9) Giupponi, G.; Buzza, D. M. A. J. Chem. Phys. 2004, 120, 10290. (10) Freire, J. J.; Rodríguez, E.; Rubio, A. M. J. Chem. Phys. 2005, 123, 154901. (11) Likos, C. N.; Schmidt, M.; Löwen, H.; Ballauff, M.; Pötschke, D.; Lindner, P. Macromolecules 2001, 34, 2914. (12) Likos, C. N.; Ballauff, M. Top. Curr. Chem. 2005, 245, 239. (13) Pötschke, D.; Ballauff, M.; Lindner, P.; Fischer, M.; Vögtle, F. Macromol. Chem. Phys. 2000, 201, 330. (14) Likos, C. N.; Rosenfeldt, S.; Dingenouts, N.; Ballauff, M.; Linder, P.; Werner, N.; Vögtle, F. J. Chem. Phys. 2002, 117, 1869. (15) Mladek, B. M.; Kahl, G.; Likos, C. N. Phys. Rev. Lett. 2008, 100, 028301. (16) Huissmann, S.; Likos, C. N.; Blaak, L. Soft Matter 2011, 7, 8419. 5387

dx.doi.org/10.1021/ma501127f | Macromolecules 2014, 47, 5379−5387