A Hybrid Monte Carlo Self-Consistent Field Model of Physical Gels of Telechelic Polymers.

We developed a hybrid Monte Carlo self-consistent field technique to model physical gels composed of ABA triblock copolymers and gain insight into the structure and interactions in such gels. The associative A blocks of the polymers are confined to small volumes called nodes, while the B block can move freely as long as it is connected to the A blocks. A Monte Carlo algorithm is used to sample the node configurations on a lattice, and Scheutjens-Fleer self-consistent field (SF-SCF) equations are used to determine the change in free energy. The advantage of this approach over more coarse grained methods is that we do not need to predefine an interaction potential between the nodes. Using this MC-SCF hybrid simulation, we determined the radial distribution functions of the nodes and structure factors and osmotic compressibilities of the gels. For a high number of polymers per node and a solvent-B Flory-Huggins interaction parameter of 0.5, phase separation is predicted. Because of limitations in the simulation volume, we did however not establish the full phase diagram. For comparison, we performed some coarse-grained MC simulations in which the nodes are modeled as single particles with pair potentials extracted from SF-SCF calculations. At intermediate concentrations, these simulations gave qualitatively similar results as the MC-SCF hybrid. However, at relatively low and high polymer volume fractions, the structure of the coarse-grained gels is significantly different because higher-order interactions between the nodes are not accounted for. Finally, we compare the predictions of the MC-SCF simulations with experimental and modeling data on telechelic polymer networks from literature.

Introduction

Here we describe a combination of the Scheutjens–Fleer self-consistent field theory with a Monte Carlo algorithm, which is used to simulate a gel network of symmetric telechelic polymers. These polymers have associative end-blocks, while the middle block is soluble. This combination leads to the formation of micelles in which the end-blocks associate in the core and the middle blocks form the corona. Such micelles are called flower-like micelles, with the core as the heart and the polymer loops as petals. At a sufficiently high concentration of micelles, the micellar cores are so close to each other that the ends of the polymers can be in different micelles, thus forming a bridge. Because the polymers can now form both loops and bridges, the number of possible conformations and thus the entropy increases. This increase in entropy gives an attractive contribution to the interaction between the micelles. If there are enough bridges to form a percolating network, a gel network is formed with the micellar cores as the nodes, as shown in Figure .

Figure 1

A network formed by triblock copolymers. The circles are the nodes, formed by the micellar cores with the associated end-blocks (indicated by gray lines), connected by the soluble middle blocks (black lines). The telechelics can form loops (both ends in the same node) or bridges (each end in a different node). An isolated flower-like micelle is shown in the lower right corner. Reproduced with permission from ref (9). Copyright PCCP Owner Societies, 2015.

Some researchers have reported that the attraction can become so strong that phase separation occurs;[1−3] others however did not observe phase separation.[4] One reason why these experiments show different outcomes is that it is difficult to synthesize these polymers. Often the middle blocks show considerable polydispersity and not all polymer ends are functionalized. The latter will increase steric repulsion between the micelles and thus prevent phase separation. In computer simulations, these problems can be avoided.

We assume that the binding energy of end-blocks to the micellar cores is so high that the concentration of free ends is negligible but still low enough that the ends can exchange between the cores. This allows the polymers to redistribute themselves over the micelles and form new bridges. This enables these gel networks to heal themselves when damaged.[5,6]

Because of these properties, telechelic polymers are applied in the paint industry to improve the rheological behavior of paints. They can also be used as a gel material for gel electrophoresis. Furthermore, they are studied as a drug carrier for slow drug release. A hydrophobic drug can be dissolved in the core of the micelles. When a gel made of telechelic polymers is placed in the body, it will slowly release individual micelles and thus the drug over time. An additional advantage in chemotherapy is that, due to the increased permeability of blood vessels in tumors, the micelles can accumulate in tumor tissue, which then receives a higher dose of the drug.[7]

There have been many experimental studies on gels made of telechelic polymers.[1−6] The number of theoretical studies and simulations is, however, limited.[8] The length of the polymers makes it time-consuming to study these networks with molecular dynamics simulations even when coarse-grained bead and spring models are used for the chains. This is because a representative fraction of all the possible states of the system has to be sampled. As the polymers are large and can easily entangle, they diffuse slowly and it therefore takes a long time to reach and sample the equilibrium structure. One could choose to use even more coarse-grained models, for example, by simulating an entire micelle as a single particle.[8] It is, however, difficult to describe the proper interaction potentials between the micelles, as the interactions are not necessarily pairwise additive. That is, the strength of the interaction between two micelles depends on the surroundings of the micelles.[9]

A solution to this problem is to employ a hybrid simulation technique which combines the benefits of particle simulations with the computational efficiency of self-consistent field calculations. Here, we combine the SF-SCF (Scheutjens–Fleer self-consistent field) method with a Monte Carlo algorithm. With the SF-SCF model, the free energy of a particular configuration of the nodes is calculated based on an average over all possible freely jointed chain conformations. This reduces the simulation time because the polymer configurations no longer need to be sampled individually. The positions of the nodes are sampled using a Monte Carlo algorithm that uses the free energy determined by the SF-SCF model to accept or reject the moves of the cores. In our model, we focus on the interactions between the nodes caused by polymers. The binding of the polymers to the nodes is done in a simplistic manner because the interactions between the nodes are not influenced by the specific mechanism through which the polymer ends bind. Hence the results can be applied to a variety of gels with polymers with associative end groups regardless of the exact binding mechanism.

The goals of this article are to demonstrate this hybrid Monte Carlo SF-SCF approach and to apply it to a system of polymer micelles in solution to gain insight in the structure of such systems. We compare the results to Monte Carlo simulations where the nodes, with polymers, have been coarse-grained to particles, with effective pair potentials as calculated in our previous article.[9] To further validate the method, the structure factors of the systems simulated by the Monte Carlo SF-SCF method are compared to experimental data found in literature.

It should be noted that this is not the first time that the SF-SCF theory is combined with a Monte Carlo algorithm. Previously, we showed some preliminary results for a charged polymer gel adsorbed on a wall, obtained with a model very similar to the one described in the present article.[10] Furthermore, Charlaganov et al.[11] used a combination of SF-SCF with Monte Carlo to study the depletion interaction of polymers near walls. They used approximate pair potentials to do a Monte Carlo simulation and subsequently used the self-consistent field theory to calculate a more accurate free energy and correct for the wrong weighing of the states. Potentially their method is more efficient, as the SF-SCF equations do not need to be solved for the rejected states. The rate at which the states of the system are visited is, however, determined by the free energy of the Monte Carlo simulation. If this free energy is not accurate, more steps are needed to reduce the noise level. This method is therefore only effective if a good approximation for the free energy can be determined. The more particles are present, the more accurate the interaction potential needs to be as the error would scale with the root of the number of particles. With our method, we do not need approximate potentials and the number of particles we could simulate is thus not limited in this way.

Method

First, we will explain the SF-SCF theory, specifically for the 3D simple cubic lattice that was used in the present study. It is based on the method used in our previous paper on the interactions between nodes with telechelic polymers, and some more details can be found there.[9] Descriptions of the SF-SCF theory for other types of lattices can be found in literature.[9,12−15] Next, we will show how we modeled the physical gel with the SF-SCF theory. Subsequently, the details of the Monte Carlo method will be described, and finally the methods for analyzing the data will be discussed.

SF-SCF Theory

With the SF-SCF method, space is divided into lattice sites, which in the present study have a simple cubic ordering. Small molecules, such as solvent molecules, are represented by a single segment that has the size of one lattice site. Larger molecules, such as the polymers considered here, are represented by multiple segments. We assume that the segments of a polymer are connected like a freely jointed chain. Because we use a simple cubic lattice, the angle between subsequent segments can only be 180°, 90°, or 0°. For 0°, the polymer is thus allowed to fold back onto itself. Segments adjacent to each other in the molecule of course still have to be next to each other on the lattice. The short-range part of the interaction between different types of segments is quantified by the Flory–Huggins parameter χ, which is half of the change in free energy when two segments are exchanged between homogeneous phases of each segment type.

It would be far too much work to generate all the ways in which the polymers can distribute themselves over the system one by one. So instead, we determine the average distribution of the polymers over the system, i.e., we try to find the volume fractions for each segment type at each lattice site. These volume fractions can also be regarded as an average over time. This is done by generating all the possible polymer conformations, which are all the possible paths of the polymer chain on the lattice. Subsequently, the polymers are distributed over these conformations according to their Boltzmann weights. Because many of the conformations are nearly identical, this saves computation time. A disadvantage is that the interactions between the segments are calculated based on the average surroundings rather than on a specific configuration of the polymers, where with a configuration we mean a particular distribution of the polymers over the conformations.

The polymers will distribute themselves over the polymer conformations according to the Boltzmann weight e– of these conformations. U is the energy in units kBT of a particular polymer conformation c, given the average surroundings of this conformation. U is the sum of the energy contributions of each segment. We call these contributions the segment potentials u(r), where X is the segment type and r is its location. These segment potentials u are calculated from the volume fractions φ of the various segment types in the neighboring lattice sites, which in turn are calculated from the segment potentials. We repeat this iterative process until we find a self-consistent solution, in other words, until the segment potentials derived from the volume fractions are the same as those that were used to calculate these volume fractions. The segment potential of a segment of type X is given by

Here, the first term ∑χ⟨φ(r)⟩ describes the average interaction energy of a segment of type X, at position r, with segments of types Y in sites adjacent to position r, χ is the Flory–Huggins parameter for the interaction between segments of type X and Y, and ⟨φ(r)⟩ is the average volume fraction of segment type Y in all neighboring lattice sites r′. The latter is given bywhere Z is the number of neighboring lattice sites. In our case, we have a simple cubic lattice and Z = 6. We do not need to consider the interaction energy between segments of the same type as the Flory–Huggins parameter χ = 0 by definition.

The second term in eq , α(r), is used to ensure that the sum of the volume fractions at each lattice site is unity. It has to increase when the sum of the volume fractions is larger than one and decrease when the sum is less than one. We chose to update α(r) with each iteration step as

The factor η = 0.3, which is small enough to prevent divergence. For the first iteration αold(r) = 0.

The volume fraction of a segment s of the polymer chain at lattice site r is given by the sum of the Boltzmann weights of all chain conformations c that pass through r with segment s, multiplied with a normalization constant C:The normalization constant C is the number of polymers n divided by the partition function q, which is the sum of the Boltzmann weights of all polymer conformations of the polymer:An efficient way to calculate q and φ is to use the propagator formalism.

The end point distribution function G(r,N + 1) is the average Boltzmann weight of all chain conformations ending with segment s = N+ 1 on lattice site r. In this study, we only allowed the polymers to start at coordinates that lie within the nodes. We therefore write the end point distribution function as G(r,N + 1|{r},0) indicating that only the conformations starting with segment 0 within {r} contribute to the end point distribution function, as we have set the Boltzmann weight of all other polymer conformations to zero. Because the position of the last segment is the same for all conformations, we can move the contribution of the last segment e– outside this summation:where X is the segment type of segment N + 1. The second part is the summation of the end point distribution functions of the chain without the last segment over all sites r′ that are adjacent to r. Because only a fraction of the conformations goes from site r′ to r, we have to multiply the propagator in r′ with . We can repeat this process until the first segment is reached. For this (starting) segment, the end point distribution function is simply e–. In this way, we can calculate the entire end point distribution function.

With these end point distribution functions, we can calculate the volume fraction of each segment s of the polymer according toAs the Boltzmann weight of segment s is in both propagators, we need to divide by e–. Because the polymers in our system are symmetric, we can save computation time by rewriting the first line of eq as the second line so that only the propagators starting with segment 0 have to be calculated.

The overall volume fraction distribution of polymers is found by summing over all the polymer segments:The distribution of the monomeric solvent S simply follows from the Boltzmann weight:When the segment potentials are normalized to zero in the pure solvent phase, which is in equilibrium with our system, C = 1.

The Helmholtz energy, which is needed for the Monte Carlo moves, is given bywhere Q is the partition function of the system. We calculate Q by using the ideal gas approximation:The first term is the contribution from the polymers, while the second term comes from the solvent. Here n is the number of solvent molecules. The single molecule partition function of the polymers q is calculated by summing the end point distribution function, over all positions r.For the solvent, the single molecule partition function is given by q = ∑e–.

In the second term of eq , we correct for the use of the Lagrange parameter. We have previously made two potentially conflicting assumptions. We have assumed that the system is incompressible, and by defining C = 1 we assume that the system is in equilibrium with a pure solvent phase. If, for example, we would place a solvophobic wall in our system, the volume fraction of the solvent would be lower near the wall than in the pure solvent. To make sure that the volume fraction next to the wall is the same as that in the pure solvent, which is required for an incompressible system, we introduced the extra potential α(r). There is of course no physical origin for this potential, and to get the correct Helmholtz energy for the given volume fractions, α(r) has to be subtracted from the Helmholtz energy.

Gel Description within SF-SCF Theory

Here and below we will express the Helmholtz energy in units kBT and measure the distances in lattice units. The polymers are represented by a chain of N = 50 segments B, which represents the middle block, and one segment A at each end, representing the end groups. We forced the end groups of these polymers to stay together, like in the micelles, by defining small volumes, called nodes, with a size of 3 by 3 by 3 lattice sites. By setting the Boltzmann weights for segments A to zero outside the nodes, the end groups are forced to stay within the nodes. The set {r} thus encompasses all lattice sites that lie within the nodes. These nodes will be moved using a Monte Carlo scheme. Because the number of nodes that we can model is limited, we use periodic boundary conditions, so there is no interface in the system.

The following values for the various parameters were used as defaults in these experiments. The Flory–Huggins parameter χ was 0.4, and the polymer volume fraction φ was 0.25. The number of nodes M was 125 with f = 5 polymers per node, thus 625 polymers in total. We investigated the effect of changing several of these parameters on the structure of the gel. The volume fraction of the polymer φ was varied from 0.5 to 0.031. The effect of the Flory–Huggins parameter was studied by doing additional calculations for χ = 0.0 and χ = 0.5. Calculations were also done for 2.5 and 10 polymers per node. The number of polymer ends in each node is not fixed but can fluctuate around the average value depending on the statistical weights of the conformations starting and ending at this node. In practice, these fluctuations in the number of polymer ends in each node are limited due to the steric hindrance between the polymers. This is similar to real systems in which the number of polymers per node can also fluctuate. It also allows for slightly different compositions of the dilute and concentrated phases when phase separation occurs. In Figure , we show a few examples of the probability density function f(Ne) of the number of end groups per node N. This distribution clearly becomes wider as the concentration increases. At high density, the steric hindrance between the polymers is less because the density around the nodes quickly drops to the bulk density and the polymers from the same node repel each other only over a short distance. It is therefore not so disadvantageous to put more than the average number of polymers on a node.

Figure 2

Probability density function f(Ne) of the number of end groups per node Ne for N = 125, χ = 0.4, and f = 5. ϕ = 0.25 (dashed line), ϕ = 0.125 (dotted line), and ϕ = 0.063 (solid line).

To see to which extent the outcome of the simulation was affected by the limited number of nodes, we also did some simulations with 8, 27, and 64 nodes. For some systems, the radial distribution function had not flattened out at a distance of half the box size. We therefore also did simulations with 512 nodes. A more detailed overview of the calculations performed can be found in the Supporting Information.

Monte Carlo Protocol

A basic Monte Carlo simulation consists of doing a Monte Carlo step, which is a trial move in the parameter space, and an acceptance rule which determines whether or not to accept the move based on the change in (free) energy. In our case, the trial moves consisted of picking a number of nodes at random and moving them by one lattice site in a random direction. A node could be selected multiple times during a single Monte Carlo step and can thus also move multiple sites. The number of nodes that are moved is adjusted during the equilibration part of the simulation, such that the acceptance ratio is about 25%. After the nodes have been moved, the distribution of the polymers is calculated again and the new Helmholtz energy Fnew is compared to the old Helmholtz energy Fold. The reason for using the Helmholtz free energy is that when a node is moved, not only the interaction energy changes but the conformational entropy of the polymers is changed as well. If the new Helmholtz energy Fnew is lower than the old Helmholtz energy Fold, the move is accepted. If it is higher, it is accepted with the probability:

At the start of the simulation, the nodes were ordered in a simple cubic ordering filling the whole cubic simulation volume. We aimed to do m = 40 000 Monte Carlo steps in each simulation. This is long enough for the system to equilibrate provided that the density remains homogeneous. To demonstrate that the system is equilibrated well within 40 000 steps, we show the SF-SCF Helmholtz energy as a function of the number of Monte Carlo steps in Figure .

Figure 3

Helmholtz energy per polymer as a function of the number of Monte Carlo steps m. χ = 0.4, ϕ = 0.25, f = 5, and M = 125.

At first sight it may seem puzzling that the Helmholtz energy increases as the system relaxes. The entropy of the nodes is, however, not included in the Helmholtz energy presented in Figure . At the start, the nodes are in a highly ordered state. By distributing themselves more randomly over the volume, the entropy of the nodes is increased. This results in a lower Helmholtz energy for the system as a whole even though the Helmholtz energy of the polymer chains has increased. In principle, the entropy of the nodes can be calculated from the radial distribution function and higher-order particle correlation functions using Green’s entropy expansion.[16] In our case, the three-particle correlation function was still rather noisy and it was therefore not possible to accurately determine the entropy of the nodes.

Coarse Grained Simulation

To show that the hybrid Monte Carlo SF-SCF method describes the system better than Monte Carlo simulations with coarse grained nodes, we performed Monte Carlo simulations with M = 125, f = 5, and χ = 0.4. In these simulations, the nodes with their polymers have been coarse-grained to a single particle. We used an effective interaction potential, calculated as described in our previous article,[9] as the interaction potential between these particles. To determine this effective pair potential, we first calculated the free energy per node for a simple cubic arrangement for different distances between the nodes. Subsequently, we calculated the effective pair potential so that it gives the correct free energy for all distances. The resulting potential is shown in Figure . The depth of the well is 0.33kBT.

Figure 4

Effective pair potential used in the Monte Carlo simulations. The distance r is measured in lattice sites.

Data Analysis

We calculated the radial distribution function of the nodes to see how much ordering there is in the system. This was done by splitting the range of possible interparticle distances in a number of subranges called bins. The width of these bins is dr. Next, we loop over all particle pairs and Monte Carlo steps m and count how many particle pairs have an interparticle distance that would fall in each bin b. ”nint” indicates that number is rounded to the nearest integer.In this equation, V is the volume in the number of lattice sites, V is the number of lattice sites that fall within the bin b(r) at radius r, and m is the number of Monte Carlo steps over which the radial distribution function is averaged. M is the number of nodes, and rn is the position of the node.

To be able to compare the results of these simulations to experiments, we also calculated a structure factor S(ξ) based on the radial distribution function using

In this equation, g(r) is the radial distribution function, r the distance, ρ the number density of the nodes, and ξ the spatial frequency.

Because of the finite size of our system, the radial distribution function does not go to exactly unity for large distances. This can, for example, be seen in Figure , where the dotted curve for M = 125 stays just above unity. The explanation is that if a particle has an excluded volume, the volume remaining for the other M – 1 particles is a bit smaller and the radial distribution function will be a little bit higher than unity far away from the particle. Similarly, if the interaction between the nodes is attractive, the concentration close to the node will be higher and far away and it will be a bit lower. In that case, the radial distribution far away will be a bit less than unity. As a result, a peak shows up around ξ = 0 in the structure factor. As the osmotic compressibility is effectively determined by extrapolating the structure factor to zero, we need a way to suppress this peak at ξ = 0.

Figure 5

Effect of the number of nodes used in the simulation on the radial distribution function. The numbers of nodes in the system are 8 (bold solid line), 27 (dashed line), 64 (solid line), and 125 (dotted line). χ = 0.4, ϕ = 0.25, and f = 5.

Recently, Dawass et. al[17] wrote an article comparing several methods for correcting some of these finite size effects. The best method according to them was the method of Ganguly and van der Veght.[18] They adjusted the radial distribution function at distance r based on the excess amount within distance r. To us this did not seem optimal, as the excess amount fluctuates considerably as a function of the distance. As a first-order approximation, the value of the entire radial distribution function will be increased due to the local excluded volume of a particle. We therefore think that a correction that is more uniform would be better at approximating the real radial distribution function. The most logical thing to do would thus be to multiply the radial distribution with a small factor such that the radial distribution function goes to exactly 1 at long distances. For small simulation volumes, the radial distribution function is however not yet entirely flat at a distance of half the box size. It is thus not so easy to determine what the right correction factor is. Ideally, we get a smooth curve near a spatial frequency ξ = 0. If we, however, get it wrong, there is a significant spike in the structure factor close to ξ = 0. It therefore seemed reasonable to choose this correction factor such that the magnitude of the second derivative near ξ = 0 is minimal, although a different derivative may work as well. To determine whether our method works, we simulated two systems, one with hard spheres and one with the effective interactions we use in the coarse grained simulation. We compared the corrected radial distribution function and compressibility of boxes with 512 particles to those of a simulation box with 13824 particles to see if our method would give a useful correction of the radial distribution function. The corrected radial distribution functions for the systems with 512 particles give Kirkwood–Buff integrals that deviate less than 15% from the value of the large system, while the uncorrected values deviated as much as 60%. For the hard sphere system with 512 particles, the value differs by about 5% from the theoretical value obtained with the K-equation of state,[19] and for the system with 13824 particles, our correction reduced the deviation from 5% to 1.5%. To our knowledge, this method has not been described in the literature and we hope to soon write a short communication in which we compare this method to other methods for correcting finite size effects. The values of the correction factors ranged from 0.987 to 1.008. With this corrected radial distribution function, we calculated the osmotic compressibility κ according to

Results and Discussion

In Figure , an example of the simulation volume is shown. The nodes are clearly visible as dark-red cubes with a slightly lighter core. Because of the steric repulsion between them, the polymers push each other away from the node and so drag their anchoring groups to the outside of the node. This results in a relatively low density within the core of the node.

Figure 6

3D view of the default system. The nodes are colored red. The polymer concentration decreases as the color goes from red via white to blue.

In Figure , we show the radial distribution function for systems with different numbers of nodes M and thus also different volumes. All other parameters have their default values. For M = 8 and M = 27, the radial distribution functions deviate significantly from the ones for M = 64 and M = 125, which are very similar. It thus seems that our default conditions using 125 nodes gives results that do not deviate too much from an infinite system, although there is still some effect of the limited size of the simulation volume. The system should not be much smaller, as the peak of the second coordination shell has barely ended at a distance equal to half the box size. With ϕ = 0.5 and f = 10, the radial distribution function shows peaks well beyond half the box size, and for this system as well as several others, we performed simulations with M = 512 nodes. This still is not optimal, but the computational costs are too high to simulate even larger systems.

The dependence of the radial distribution function on the overall polymer volume fraction is shown in Figure . As the polymer volume fraction is increased from ϕ = 0.125 to ϕ = 0.5, the peak of the radial distribution function shifts inward. Hence, at high concentrations, the polymers are pressed into each other as there is not enough space to place all nodes at their optimal distances. As the volume fraction is reduced, the distances between the nodes increase until the optimal distance is reached at a volume fraction of about ϕ = 0.125. There is no strong ordering in the sample, and only two relatively weak coordination shells are visible in the radial distribution function. The system is thus expected to behave in a liquid-like manner on time scales that are long compared to the relaxation time of an individual bridge.

Figure 7

Effect of the polymer volume fraction on the radial distribution function. The tested volume fractions are ϕ = 0.50 (bold solid line), ϕ = 0.25 (bold dashed line), ϕ = 0.125 (bold dotted line), ϕ = 0.063 (solid line), ϕ = 0.031 (dashed line), and ϕ = 0.0078 (dotted line). χ = 0.4, M = 125, f = 5.

At first sight, it may be surprising that the level of ordering of the nodes does not increase with increasing node concentration. One would expect that due to the strong steric repulsion the nodes would order themselves. Similar to polymeric solutions, however, the environment starts to look more like a polymer melt, as the polymer concentration is increased. The polymers are therefore distributed more homogeneously over the volume. As a result, the steric hindrance experienced by the nodes will depend less on their position and the system can thus remain unordered.

As the polymer volume fraction is decreased from ϕ = 0.125, the peak of the first coordination shell rises, suggesting that the strength of the attraction increases. This is in line with our earlier finding that the interaction between two nodes depends on their surroundings.[9] As the system is diluted, the number of neighboring nodes decreases and the attraction with the remaining neighbors increases. For dilute systems, the binding energy can be estimated by taking the logarithm of the peak height of the radial distribution function. In this case, the height is about 2.7 for ϕ = 0.0078, which corresponds to a binding energy of roughly 1 kBT. This binding energy and the position of the peak are the same as we found before.[9]

Let us next consider the effect of solvent quality. In Figure , radial distribution functions are shown for different values of χ. At a volume fraction of ϕ = 0.50, shown in Figure a, there is practically no difference between the different solvent qualities. At such a high polymer volume fraction, the swelling of the polymer corona does not significantly decrease the number of unfavorable interactions between the polymer segments because they would swell into the corona of the next micelle. The size of the micelles in a good solvent is therefore the same as that in a theta solvent, and the radial distribution function is therefore also practically the same.

Figure 8

Effect of the solvent quality χ on the radial distribution function for different volume fractions: (a) ϕ = 0.50, (b) ϕ = 0.25, (c) ϕ = 0.125, and (d) ϕ = 0.0625. χ = 0.5 (bold solid line), χ = 0.4 (bold dashed line), and χ = 0.0 (dotted line). M = 125, f = 5. In parts c and d, part of the radial distribution lies beyond half the box size. The values in this range are displayed to show that phase separation occurs, although they are probably still be affected by the limited box size.

This is illustrated in Figure , where the interaction potential ΔF between two isolated nodes is plotted for different background polymer concentrations. At a background polymer volume fraction of ϕ = 0.5, the curves for χ = 0.0 and χ = 0.5 are practically the same.

Figure 9

Effect of background polymers on the interaction potential between two nodes: χ = 0.0 black and χ = 0.5 gray. ϕb = 0 (dotted line), ϕb = 0.1 (bold dashed line), and ϕb = 0.5 (bold solid line). f = 5.

As the polymer volume fraction is decreased to ϕ = 0.25, the radial distribution functions for the different solvent qualities start to differ. The peaks of the radial distribution functions shift outward, most strongly for the good solvent. For the theta solvent, the radial distribution function is otherwise similar to that at ϕ = 0.50. For a good solvent, the height of the peaks increases, as the steric repulsion is strongest in a good solvent and it thus gives the most ordered structure.

When the volume fraction is lowered further to ϕ = 0.125, we observe that for χ = 0.5 the radial distribution function no longer goes to unity at large distances. This is most likely because phase separation occurs: the cross section of the gel in Figure clearly shows a dilute and a concentrated region.

Figure 10

A cross section of the gel at χ = 0.5, ϕ = 0.125, f = 5. The darker the color, the higher the polymer density. The maximum polymer volume fraction is about 0.75. The nodes have clearly clumped together forming a dense region, which coexists with a dilute region with just a few micelles.

In addition, the first peak for χ = 0.5 is higher than the peak for χ = 0.0. For χ = 0.5, the interactions are now attractive and they become stronger as the gel becomes more dilute, while for χ = 0.0 there is still a net repulsion between the nodes which decreases as the gel becomes more dilute.

Finally, in Figure d, the polymer concentration has been lowered to ϕ = 0.0625 and now the radial distribution function does go to 1 for χ = 0.5. This, however, does not mean that the system is already below the lower binodal. At the start of the simulation, the nodes are distributed homogeneously over the volume. They will initially clump together in small clusters. These clusters, however, diffuse much slower than individual nodes. It will thus take a long time before all the clusters and nodes have aggregated by diffusion and Ostwald ripening and formed a dense phase. The simulation was therefore too short to fully equilibrate the system. The radial distribution function does show a slight dip at a distance of about 45 lattice sites, which is also visible for χ = 0.4. As the individual nodes and clusters diffuse around, they stick to other clusters. The concentration of micelles and other clusters near this cluster therefore decreases, leading to a zone with a relatively low concentration. This process may not only happen for complete phase separation but also in the case that the clusters have not yet reached their equilibrium size distribution. This is illustrated by the change in the radial distribution functions as the number of Monte Carlo steps is increased. The more Monte Carlo steps have been taken the further out the dip lies and the deeper it becomes. The system is thus not equilibrated within the simulated number of Monte Carlo steps. At the end of the simulation, there is also a clear void visible within the gel.

It is possible to improve the rate at which the system equilibrates by occasionally making Monte Carlo moves that displace micelles over large distances. We, however, intended to study the homogeneous phases of these micellar solutions and therefore did not implement such large Monte Carlo moves.

At a volume fraction of ϕ = 0.0625, the distance between the nodes is so large that, for all χ, the interactions are no longer repulsive at the average intermicelle distance. The peak for χ = 0.4 is therefore higher than that at χ = 0.0 because the height of the peaks is now determined by the strength of the attraction between the micelles.

Now we turn to the effect of the number of polymers per node f, as shown in Figure . For f = 2.5, the radial distribution function has just one peak just like a gas. In contrast, there are many peaks visible for f = 10. For the highest concentration ϕ = 0.5, these peaks occur beyond half the box size. It is therefore likely that in this case the radial distribution function is still influenced by the size of the simulation volume. A striking difference between f = 10 and the lower functionalities is that the height of the peaks increases as the concentration is increased from ϕ = 0.25 to ϕ = 0.50. This suggests that as the number of polymers increases, the micelles will behave more like hard particles and crystallization may be possible for nodes with even more polymers.

Figure 11

Effect of the number of polymers per node f on the radial distribution function for different volume fractions: (a) ϕ = 0.50, (b) ϕ = 0.25, (c) ϕ = 0.125, and (d) ϕ = 0.0625. f = 2.5 (bold solid line), f = 5 (bold dashed line), and f = 10 (bold dotted line). χ = 0.4, M = 125, except for f = 5 with ϕ = 0.5 and f = 10 with ϕ = 0.5 or ϕ = 0.25, where M = 512.

Figure d shows that, for f = 10, the radial distribution function drops a bit below unity at large distances, although the deviation is not as large as in Figure c. There are some interconnected cavities visible within the gel. It is therefore possible that this gel will also undergo phase separation even though this is not yet clearly visible. The number of Monte Carlo steps taken is relatively small, and the gel may not have had enough “time” to phase separate.

Now that we have discussed the radial distribution functions for different parameters, we can compare them with the radial distribution functions calculated with Monte Carlo simulations in which we coarse-grained the nodes as single particles. In Figure , radial distribution functions from the MC-SCF simulations and the Monte Carlo simulations with effective pair potentials are shown.

Figure 12

Comparison of the radial distribution functions for the MC-SCF (bold solid line) and the normal MC simulation with effective pair potentials (dotted line). The different volume fractions of the polymer are (a) ϕ = 0.50, (b) ϕ = 0.125, and (c) ϕ = 0.031. In all simulations, f = 5; in (a), 512 particles were used, and in (b) and (c), 125.

At high densities (see Figure a), the MC simulation with effective pair potential gives much sharper peaks than the MC-SCF model. This is probably caused by an overestimation of the repulsive force between the particles. When two nodes approach each other closely, the polymers can move out of the way if there are no other particles nearby. However, if the nodes have many close-by neighbors, the polymers can not move out of the way and the repulsion is thus stronger. The MC-SCF model can distinguish between these cases. A pairwise interaction, however, cannot, and instead an assumption has to be made about the surroundings of the nodes. In the way we determined the effective pair potential, it is assumed that the other nodes are at the same distance from the interacting nodes as the interacting nodes are from each other. When a particle is closer than the typical distance between a particle and its nearest neighbors, the average distance to the other nodes is underestimated and the repulsive force is too strong. Because of this, the nodes cannot approach each other as closely as in the MC-SCF model and therefore appear as harder particles, resulting in the sharper peaks.

At low concentrations, as seen in Figure c, the opposite problem arises. Here the peak of the first coordination shell is much higher for the MC-SCF model. Not only does the effective pair model overestimate the repulsion between the particles, it also underestimates the attraction. When a node already has many neighbors, adding another one increases the number of polymer conformations relatively little compared to the total number of potential conformations. If instead a node has no neighbors, the relative increase in the number of polymer conformations is much larger. The change in free energy when a neighbor is added will therefore be larger when a node has fewer neighbors. The attraction at low concentration will therefore be stronger. With the effective pair potential, it is assumed that there are neighboring nodes at the same distance as the interacting nodes. This results in an underestimation of the attraction at low concentration.

At intermediate concentrations (Figure b), the effective pair interaction gives roughly the same radial distribution function as the MC-SCF model, although the repulsion between the micelles is still overestimated at short ranges.

On the basis of these calculations, it is clear that a Monte Carlo simulation with a single pair potential does not correctly describe the behavior of the nodes at a wide range of concentrations, although some improvement should be possible as the short-range repulsion appears too strong at all concentrations. An option would be to use a custom potential for each density. For systems in which the density remains homogeneous, this would be an improvement. If the density is, however, not homogeneous, the result would be even worse than with the effective pair potential we have used here.

To validate our MC-SCF simulations, we need to make predictions that can be compared to experimentally obtained results. We therefore determined the structure factor and the osmotic compressibility κ according to eqs and 17.

The structure factors are shown in Figure and the compressibility is plotted in Figure . Close to spatial frequency ξ = 0, the uncertainty in the structure factor is relatively large. As the structure factor near ξ = 0 is closely related to the compressibility, the accuracy with which the compressibility can be calculated is also limited.

Figure 13

Structure factors calculated from the radial distribution functions for: (a) χ = 0f = 5, (b) χ = 0.4, f = 5, (c) χ = 0.5, f = 5, and (d) χ = 0.4, f = 10. ϕ = 0.50 (bold solid line), ϕ = 0.25 (bold long-dashed line), ϕ = 0.125 (bold dashed line), ϕ = 0.06 (bold short-dashed line), and ϕ = 0.03 (dashed line). The gray areas under the graphs indicate the 99% confidence interval.

Figure 14

Osmotic compressibility relative to that of an ideal gas. (a) For different numbers of polymers per node: f = 2.5 (bold solid line), f = 5 (bold dashed line), and f = 10 (bold dotted line). (b) As a function of χ: χ = 0 (bold solid line) and χ = 0.4 (bold dashed line). The error bars indicate the 99% confidence interval.

In two of the presented cases, the structure factors are negative at ξ = 0. For systems in equilibrium, this is physically unrealistic and it most likely results from the limited size of our simulation volume. In several other cases, increasing the number of nodes from 125 to 512 caused the structure factor to become positive. This would probably also be the case for these two systems if we could run the simulations with more Monte Carlo steps and nodes.

The effect of the number of polymers per node f on the osmotic compressibility (eq ) is shown in Figure a. The values shown are relative to the compressibility of an ideal gas with a particle concentration that is the same as the concentration of nodes in our simulations. At high polymer volume fractions, the steric repulsion of the polymer coronas prevents the nodes from coming close to each other and the osmotic compressibility is therefore much smaller than that of an ideal gas. At low concentration, the attraction causes the nodes to form clusters and the osmotic compressibility is therefore higher than that of an ideal gas because the number of freely moving particles is reduced. The relative compressibility gives a lower limit for the number of nodes that form a cluster. For χ = 0.4, f = 5, and ϕ = 0.03, the average cluster size should be, for example, at least 3.

For ϕ = 0.5, the nodes with the fewest polymers per node seem to have the highest relative compressibility. The confidence intervals, however, still overlap so the difference is not significant. If we had instead looked at the real compressibility, the order would probably be reversed. The higher the number of polymers per node, the higher the concentration of polymers close to the nodes, while the polymer concentration halfway between the micelles is lower and thus the steric repulsion will be lower as well.

As the concentration is lowered, the order quickly changes. Because the nodes with more polymers have more attraction between them and thus form larger clusters, these systems are more compressible at low volume fractions.

In Figure b, the compressibility as a function of the Flory–Huggins parameter χ is shown. As expected, the system with χ = 0 is the least compressible; as the corona swells the most, the steric repulsion is the strongest for this case. Above, we concluded, on the basis of the radial distribution function, that phase separation occurs for the combination of f = 5 and χ = 0.5. According to theory, the compressibility should therefore go to infinity. Our system, however, has a limited number of particles and therefore the value the compressibility can reach is limited. Furthermore, we used the entire radial distribution function to calculate the compressibilities. This, however, includes the part of the radial distribution function far away from the particle where it is below unity. This lowers the calculated value of the compressibility even further. The values obtained for χ = 0.5 and ϕ = 0.13–0.03 are thus incorrect and therefore not shown in Figure b.

One of the experimental studies in literature to which we can compare our results is by Filali et al.[1] They investigated a system of swollen surfactant micelles to which they added PEO polymers with hydrophobic end groups. Under the conditions used, the Flory–Huggins parameter for PEO is between χ = 0.4 and χ = 0.5.[20,21] Although the polymers had about 120 Kuhn segments and were thus longer than the polymers we simulated, our results should show a fairly good match, as the effective pair potential is almost identical for 50 and 100 segments when the distance from the core is rescaled.[9] In addition, the core of the micelles is larger than our nodes. However, compared to the volume of the coronas, the cores are still relatively small. The experimental results should therefore be in between our results for χ = 0.4 and χ = 0.5.

Filali et al.[1] observed phase separation for f ⩾ 6, which is not much higher than the f = 5 for which we observed phase separation with χ = 0.5. In this respect, their results fit nicely with our findings.

The authors did not report the structure factor separately but did show the total scattering intensity. As the form factor goes to unity for small ξ, we should be able to make a qualitative comparison between our structure factor and their scattering intensities at small ξ. To get realistic values for our spatial frequency, we need to choose a value for the lattice size. We chose a value of 7.4 Å because this coincides with the Kuhn length of PEG.[22]

On the basis of the Daoud Cotton model,[23] the system of Filali et al. with an oil droplet volume fraction of 7% should best match our simulations with a polymer volume fraction of ϕ = 0.03. When going from large ξ to small ξ, there is a dip after the peak, indicating the average distance between the nearest neighbors followed by a steep increase in both cases. These features are less pronounced than in our simulations. This is probaly because we used f = 5, while the experimental system had on average four polymers per micelle. For higher concentrations, the structure factors do not match because the surfactant micelles in the experimental system are charged and repel each other at these concentrations.

François et al. published several experimental articles on telechelic polymers with PEO middle blocks.[24,25] In contrast to us, they found cubic phases at high concentrations using X-ray scattering.[24] Because the number of polymers per micelle was not reported, a one-on-one comparison with our simulations is difficult. Probably the number of polymers per micelle is higher than in our simulations. This may explain why they observed a cubic phase and is corroborated by the fact that for longer middle blocks, for which the number of polymers per micelle is lower, crystallization was not found.

Another factor that may have contributed is that at least 10% of their polymers had only one functionalized end, which increases the repulsion between the micelles. They also observed phase separation for systems with relatively short middle blocks (PEO Mw ⩽6000 g/mol) but not for polymers with long middle blocks (Mw ⩾10000 g/mol). As explained before, systems with longer polymers have fewer polymers per micelle and therefore less entropic attraction due to bridge formation.

Sprakel et al. studied the rheological and phase behavior of solutions of the same type of telechelic polymers, both experimentally[4] and with computer simulations.[8,26] In contrast to our simulation, they did not observe phase separation in their experimental study.[4] Although the number of polymers per micelle was not reported, it was probably larger than 6, which was estimated by Filali et al. to be the lower boundary for phase separation. A possible explanation for not observing phase separation is that about 10% of the polymers had only one associative end group. This increases the steric repulsion between the micelles, and the net attraction, which causes the phase separation, is thus reduced.

In a second paper, Sprakel et al.[26] addressed the phase behavior of the system with a SF-SCF model in which the micelle is modeled in a 1D spherically symmetric system with a reflecting boundary condition. The number of polymers per micelle was not fixed, but instead the grand potential was optimized to determine f. Phase separation was predicted for all the combinations of middle block and end-block lengths they used in their study. The minimum number of polymers per micelle they found was about eight, but their polymers were much longer than those in our study. As they used χ = 0.5, all their systems lie above the line from f = 10 with χ = 0.4 to f = 5 with χ = 0.5. Their predictions are thus in line with what we found here.

In a third study,[8] they coarse-grained the micelles to single particles. In this case, no phase separation was found. They do not mention the Flory–Huggins parameter, but they wanted to reproduce the above-described experimental system[4] and the value of χ should thus be between 0.4 and 0.5. Given that the simulated micelles have f = 25 polymers each, phase separation would be expected based on our results. However, the interaction potentials in the coarse-grained model do not take into account that the attraction will increase with a decreasing number of neighbors.[9] Moreover, a relatively small well depth of 0.38 kBT was used, comparable to the well depth we found for f = 5 and χ = 0.4 (about 0.34 kBT). Considering that the well depth for an isolated pair of micelles roughly scales with f0.5,[9] the well depth expected for this system with f = 25 would be about twice as large. These combined factors explain why they did not find phase separation from the coarse-grained modeling.

Conclusion and Outlook

We successfully combined a Monte Carlo algorithm with the Scheutjens–Fleer self-consistent field theory. With it, we were able to calculate the radial distribution function, structure factor, and compressibilities for solutions and gels of ABA triblock copolymers with varying properties over a range of densities. For f ≤ 5 polymers per node, we found, somewhat counterintuitively, that as the polymer volume fraction ϕ increases from ϕ = 0.25 to ϕ = 0.5 the amount of ordering in the system is decreased. We argue that this is because at high volume fractions the background concentration of the polymers of the other nodes becomes more homogeneous. The amount of steric repulsion therefore depends less on the position of the node. We further discovered that for χ = 0.5 and f ≥ 5, phase separation occurs. We were, however, not able to determine the compositions of the coexisting phases as the number of simulated particles was small and there should thus be a considerable effect due to the interface. Simulating such a large volume that the effects of an interface would be negligible would take far too much computation time. To avoid the effect of the interface, the Gibbs ensemble[27] can be used. Because we used a lattice, we can, however, not change the volume by arbitrary small steps but only by one lattice layer at a time. The larger the simulation volume, the larger the change in volume and thus in free energy will be. The chance that an exchange in volume would be accepted would therefore become smaller for larger and larger systems. This limits the system size we can use in combination with the Gibbs ensemble. Instead, it may be possible to “simulate” a Gibbs ensemble by simulating two volumes which would be representative for larger volumes of the simulated Gibbs ensemble. By moving particles in and out of the simulated volumes, the density could be adjusted to that of the simulated volumes of the Gibbs ensemble. To our knowledge, such an approach has not been described in literature yet. Another approach would be to coarse-grain the micelles while maintaining the dependence of the interacting potential on the surroundings of the interacting micelles. This method would also enable us to study dynamic properties of the system.

The limited system size may have affected some of our simulations at high polymer volume fractions where the radial distribution function had not completely flattened out by half the box size. The next generation of GPUs, however, promises to have 10 times more computation power as those we used. This allows larger simulation volumes and more Monte Carlo steps, making a Monte Carlo SCF hybrid model a feasible tool for future studies. By comparing the results of coarse-grained Monte Carlo simulations with those of the hybrid MC-SCF model, we have shown the shortcomings of using only one pair potential to describe the interactions between the nodes. The MC-SCF hybrid method is therefore a useful tool to model systems of flower-like micelles and telechelic networks over a wide range of concentrations.