Dissipative particle dynamics simulations of tri-block co-polymer and water: Phase diagram validation and microstructure identification

In this study, the phase diagram of Pluronic L64 and water is simulated via dissipative particle dynamics (DPD). The peculiar structures that form when the concentration varies from dilute to dense (i.e., spherical and rod-like micelles, hexagonal and lamellar phases, as well as reverse micelles) are recognized, and predictions are found to be in good agreement with experiments. A novel clustering algorithm is used to identify the structures formed, characterize them in terms of radius of gyration and aggregation number and cluster mass distributions. Non-equilibrium simulations are also performed, in order to predict how structures are affected by shear, both via qualitative and quantitative analyses. Despite the well-known scaling problem that results in unrealistic shear rates in real units, results show that non-Newtonian behaviors can be predicted by DPD and associated with variations of the observed microstructures.


I. INTRODUCTION
Structured fluids are colloidal dispersions typically obtained by mixing an organic phase with an aqueous one with the help of surfactant molecules. Within the fluid, microor nano-phases characterized by well-defined microstructures can be obtained, by varying the nature and the concentration of the components, as well as the mixing rate. The simplest example of structured fluids, which can produce a rich variety of microstructures, is the mixture formed by amphiphilic surfactants and water. As well-known, surfactant molecules self-assemble together in water, with the hydrophilic part creating a shell around the hydrophobic core and forming different microstructures, ranging from spherical and cylindrical micelles to hexagonal and lamellar structures. Phase diagrams are built to forecast when a specific microstructure is formed, and their derivation (from scattering and rheological experiments) is a quite standard procedure in this research area. However, when a structured fluid is deformed, the fate of the involved microstructures is less explored and this work aims at addressing this issue by using a computational model.
The self-assembly of large surfactant molecules, as well as the effect of shear on the observed microstructures, takes place on time scales which is commonly not accessible by traditional All-Atom Molecular Dynamics (AAMD). Therefore, mesoscopic models, such as Dissipative Particle Dynamics (DPD), can be used to investigate wider time scale ranges and highlight peculiar behaviors of such fluids, even though the molecular resolution of the model is reduced. 1 DPD describes the interaction between beads, representative of clusters of atoms and molecules, using a bead-spring model, repulsive soft potentials combined with stochastic and dissipative forces. Prhashanna et al., 2 Cheng et al., 3 Zhen et al., 5 Cao et al., 4 and Li et al. 6 already discussed and validated the reliability of DPD to predict the formation of micelles, at thermodynamic equilibrium, and microstructures in systems composed of water and surfactants. Also, complete phase diagrams can be obtained for ternary compounds as described by Wang et al., Son et al., and Yuan et al. [7][8][9] When a mechanical perturbation is applied (like in a mixing tank), shear stresses induce deformation of the microstructure that can lead to phase transitions. [10][11][12] One of the most popular structured fluids investigated in the literature, both via experiments and simulations, [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] is composed of water and the triblock co-polymer of polyethylene oxide (PEO)-polypropylene oxide (PPO)-polyethylene oxide (PEO), under the commercial name Pluronic®, by BASF. This nonionic co-polymer can be manufactured by varying the number of EO and PO monomers such that the length of the hydrophilic (PEO) and hydrophobic (PPO) blocks can be varied to tune its amphiphilic properties, hence its phase diagram together with the temperature. Spherical and cylindrical micelles are present at low concentrations, soft-gels are formed at slightly higher concentrations, and oriented lamellar structures are observed at very high concentrations. 37,38 Due to their incredibly wide range of applications and their relative simple structure, different Pluronic copolymers, named with different labels based on the length of the repeating units, such as P84, L64, F127, and P123, have been investigated from the computational point of view using DPD under both equilibrium and non-equilibrium simulations. 2,[38][39][40] Shear effects on microstructures can be proved by the change in the apparent viscosity shown by water-surfactant systems at different shear rates as described by Newby et al. 13 and Youssry et al., 14 where the rheological properties of Poloxamer 407 in aqueous solutions were investigated. Phase transitions, changes in orientation, deformation, and coalescence of micelles are examples of the phenomena involved when these fluids are subjected to shear. 15 The transition from spherical to cylindrical or worm/rod-like micelles is also a very common event. These elongated structures can create a structured network or align themselves according to the direction of the flow such that the overall behavior of the fluid is non-Newtonian. At high concentrations, the effect of the shear on complex systems plays an even more relevant role. For example, Gentile et al. demonstrated that lamellar phases could rearrange their shape producing multilamellar vesicles, with a resulting non-Newtonian behavior. 16,17 This work aims at predicting a full phase diagram for the system described, highlighting different morphologies and quantitatively assessing their shape and number, via a clustering algorithm (see the supplementary material). Also, we investigated the effect that an applied shear has on the equilibrium microstructures and its relationship with the changes in the resulting rheological behavior. This was done here by simulating a specific Pluronic co-polymer known as L64 [i.e., (EO) 13 (PO) 30 (EO) 13 ], using DPD, for which a set of parameters, optimized against the experimental equilibrium phase diagram, are available. 41 The shear flow is simulated by means of non-equilibrium DPD simulations, and its effect on the observed morphologies is quantified via a cluster analysis that enables counting of the polymer structures and identification of their geometries. Finally, the apparent viscosity of the structured fluid is calculated and its variation with the liquid microstructures is explained.

II. MODEL DESCRIPTION
The mesoscopic technique named Dissipative Particle Dynamics (DPD) was first introduced by Hoogerbrugge and Koelman as an alternative to lattice-gas automata schemes. 42 The original model was corrected and improved by Groot, Warren, Español, and Pagonabarraga. 43-45 Among its wide range of possible applications, DPD can be used to simulate shear effects on complex fluids due to the peculiarity of preserving hydrodynamic interactions. 2,46,47 In a DPD simulation, interacting beads are representative of clusters of atoms or molecules. The interaction between beads can be described with Langevin dynamics where f i is given by the force acting on each DPD particle i, f i , is the sum of a conservative, dissipative, and stochastic term. The conservative force can be described as follows: where a ij represents the conservative soft potential parameter, r ij = r i − r j is the relative distance between two beads i and j, r ij = r ij , andr ij = r ij |rij | , and r c is the cut-off radius, a characteristic length. The dissipative force is described as follows: where γ represents the dissipative coefficient acting as an artificial drag on the beads, w D is a weight function that defines the maximum range of application of the force, and v ij is the relative velocity between two beads i and j. Its dependence on the velocity of the beads allows DPD to act as a thermostat in regulating the temperature of the system. 1 The stochastic force can be described as follows: where σ is the stochastic coefficient, w R is again a weight function, ζ ij is a random fluctuating variable with zero mean and unitary variance, and ∆t is the simulation time step. The weight functions and the stochastic and dissipative coefficients are connected by the fluctuation-dissipation theorem 1 as follows: where k B is the Boltzmann constant and T is the temperature of the system. This last equation clarifies that the choice of one of the interaction parameters implies that the other is already defined. Bonded interactions are needed to maintain the topology of the polymer chain. Two types of bonded potentials have been investigated in this study: harmonic and finite-extensible nonlinear-elastic (FENE) potentials. The harmonic potential is described as follows: where K Harm is the harmonic constant and r e is the equilibrium distance between two connected beads, while the FENE potential is where K FENE is the FENE bond constant, r e is the equilibrium distance, and r is the distance between two beads. The viscosity of a DPD system has been calculated using the non-equilibrium method known as Lees-Edwards boundary conditions (LEBCs), 48 where different values of shear stress can be obtained through the application of different velocities on the beads that are close to the boundaries (top and bottom) of the simulation box. The maximum value of the velocity at the top of the box is equal to ½γl, whereγ is the shear value imposed on the system and l is the length of the box. If the conservation of momentum is respected, a linear velocity profile, across the simulation box, is obtained. The magnitude of the shear stress should lead to velocities that are larger than the thermal velocity of the beads, leading to meaningfully observable shear flows in computational studies. Also, periodic boundary conditions were used to ensure that the total number of beads and their behavior were consistent with the streaming effect. The viscosity can then be obtained via the following relationship: whereγ is the imposed shear rate, while P xy is the nonzero, with xy being the non-diagonal component of the stress tensor. LEBC allows us to obtain rheograms, linking any changes in the fluid morphology, due to the application of the shear, to changes in its viscosity. It must be however highlighted that, for a fixed value of the dissipative coefficient, γ, one drawback of the method is that only a limited range of shear rates is applicable without obtaining anomalies in the viscosity, as it is possible to observe in the supplementary material.
A final remark regards the conversion between DPD and real units. In order to compare the results of DPD simulations with experiments, it is necessary to define a conversion benchmark, such as a set of values representing physical quantities. Having in mind that DPD beads are representative of different clusters of atoms/molecules with the same size and weight, three variables can be used to define one possible conversion set (i.e., a length, a mass, and a kinetic energy). The length, defined as a cut-off radius, represents the maximum level of interaction between DPD beads, the mass represents the number of particles clustered into one bead, and the kinetic energy is an indicator of the thermal velocity of the beads. DPD simulations are performed using these parameters normalized to unity. If this set is fixed, all the remaining parameters can be obtained by their combination. Although this conversion set of parameters is producing consistent results in evaluating equilibrium properties, the same does not apply in non-equilibrium conditions. The conversion of DPD values into real physical units, according to the equilibrium conversion set, could produce unrealistic values for non-equilibrium quantities, for example, the actual shear rate applied (similar to what happens in non-equilibrium AAMD). It is also necessary to say that, given this set of parameters and types of interactions, chain-crossing is allowed. This could also result in deviations from real, physical quantities. Despite all these limitations we think, the present analysis is useful and can generate interesting results. In particular, this work is focused on the calculation of physical quantities, but on morphological changes that appear in the microstructures of structured fluids, when they are affected by shear stresses and on how they can be assessed from a quantitative point of view through a population analysis.

III. SIMULATION DETAILS
The simulated systems, composed of DPD beads, represent water molecules and Pluronic L64 chains. Two sets of simulations have been performed, in order to assess both equilibrium and non-equilibrium properties. The computational code used for this purpose was LAMMPS 49 (Large-scale Atomic/Molecular Massively Parallel Simulator), and graphical outputs were produced using Visual Molecular Dynamics (VMD). 50 The DPD parameters used to simulate the polymer chains were obtained from the literature. 2 The level of coarse-graining adopted to describe the Pluronic L64 chains is 4.3 for the EO repeated units and 3.3 for the PO repeated units. This means that one coarse-grained bead of EO contains 4.3 atomistic EO monomers and the same conversion procedure applies for PO. Using this set of parameters, Pluronic L64 chains are composed of 15 beads and simulated as A 3 B 9 A 3 DPD chains, where A is the coarse-grained bead for the EO unit and B is the coarse-grained bead for the PO one. Simulations of different concentrations of Pluronic L64 in water were performed by varying the number of beads of the two components, keeping the total number of beads in each box fixed (e.g., for a system composed of 81 000 beads, if 50% is composed of water, 40 500 spherical beads are water-type). Bonded and non-bonded interactions between beads are accounted for in the DPD model. The former were described using both harmonic and FENE potentials, while the latter are reported in Table I. All values are reported in DPD units.
The dissipative parameter γ was set equal to 4.5 (in DPD units) for all the species, while the stochastic parameter σ was set equal to 3, according to the fluctuation-dissipation theorem 1 when the value of k B T is equal to 1. Simulation boxes of different lengths were tested, from 20 × 20 × 20 cutoff radii to 40 × 40 × 40 cut-off radii. The shape of the box was not changed; i.e., only cubic shapes were used, in order to reproduce a bulk system surrounded by its replicas where properties are not affected by the direction of the application of external perturbations. The simulation box of 30 × 30 × 30 cut-off radii was found to be a reasonable compromise between reduction of simulation box artifacts and acceptable simulation times. The initial configuration of the system is prepared by random positioning of water beads and Pluronic L64 chains. The number density (i.e., the number of beads per unit volume) was set equal to 3 DPD units, 1 meaning that the total number of beads was 81 000.
The cut-off radius for non-bonded interactions was set equal to 1 with a time step of 0.01 DPD units. Equilibrium only simulations were carried out for 2 × 10 6 time steps, while non-equilibrium simulations were carried out for 3 × 10 5 followed by 5 × 10 5 time steps applying different shear rates in different simulations. The range of concentrations spans from 5% to 95% in weight percentage (w/w) of Pluronic L64, and the range of non-dimensional DPD shear rates varies from 0.005 to 2. The DPD energy was set equal to 1, and its value was recorded every 500 time steps. Moreover, the DPD property of preserving hydrodynamic interactions ensures that momentum is conserved across the box. The modified velocity Verlet algorithm, according to the scheme proposed by Groot and Warren, 1 was used as an integration scheme. During the overall simulation time, the energy was stable at 1.0 ± 0.01 k B T. A cluster analysis algorithm was written in Python 51 (Python Software Foundation, https://www.python.org) in order to quantify modifications in the microstructure when shear is applied. The code uses the Density-based spatial clustering of applications with noise (DBSCAN) 52 and it has been modified and adapted to the output file generated by LAMMPS such that clusters can be detected, and also gyration radius and cluster mass distribution (CMD) plotted in time. A complete description of the code can be found in the supplementary material (code can be downloaded at: https://github.com/hermessc/clusteringAlgo). The python code takes coordinates in input every 500 time steps and a density-based algorithm captures the number of clusters/structures of the central hydrophobic part of the Pluronic L64 chain, while the other beads (i.e., water and PEO) are ignored. In such a way, it is possible to find specific patterns in the microstructures, such as formation of spherical micelles, soft-gel, hexagonal phase, and lamellar structures. In equilibrium and non-equilibrium simulations, clusters are recognized according to a cutoff distance between closer beads. Closer neighbors are assigned to one single cluster, and beads that are outside of the cutoff range are considered as belonging to other clusters. The expected total number of clusters is not known a priori. The number of clusters was therefore monitored against the simulation time (for over 2 × 10 5 time steps) in both the equilibrium and non-equilibrium stages. Also, the cluster size, quantified by the aggregation number, N, namely, the number of Pluronic L64 chains in one cluster, was monitored, and snapshots of the configuration of the system were stored. These data were then used to identify the cluster mass distribution (CMD) indicating how the population of clusters in the simulation box is distributed over the aggregation number N.
In the lower range of concentrations, the data collected during the cluster analysis were employed to calculate the micelles gyration radius and were used to determine their sphericity. It must be highlighted that, for complex structures, there is a correlation between the aggregation number and the gyration radius N = CR d g , (12) where N is the already introduced aggregation number, C is a constant, R g is the gyration radius of the micelle, and d is a scaling exponent that tends to three in the case of spherical structures and tends to two in the case of cylindrical or worm-like structures. Plotting the aggregation number versus the radius of gyration (or vice versa) in a log-log scale allows us to identify the value of the exponent d.
In non-equilibrium simulations, each system was initialized with a linear velocity profile, with the maximum desired velocity at the top of the box and zero velocity at the bottom (see, for example, Fig. S2 of the supplementary material).
Shear was only applied on the xz plane, meaning that only P xy , one of the three non-diagonal components of the stress tensor, was not null. In LAMMPS, LEBC is not directly implemented and the fix deform function can be used for the purpose. This fix ensures that the box is deformed at a constant shear rate and together, with the remap option, positions and velocities of the particles crossing the boundaries are stored and modified according to LEBC (i.e., an extra path is added when a particle crosses a periodic boundary due to the different velocities on the top and bottom of the box). If the applied deformation exceeds a limit value such that the simulation box becomes too skewed, the domain is tilted, and particles remapped into an equivalent configuration. A linear velocity profile, hence a constant shear rate, can be obtained in this way. By using this fix, the velocity of the bottom slab is equal to zero, while the velocity on the top slab becomesγl, whereγ the DPD shear rate value and l is the length of the box. In this operative regime (i.e., from 0.005 to 2 DPD shear rate), a linear velocity profile was obtained for both systems containing only water and a mixture of water and Pluronic L64.
In order to ensure the validity of the results in the operative range, two sets of tests were performed for the upper and lower limits of the shear range. To set the upper limit, we observed the behavior of water viscosity, which needs to be consistent with its Newtonian nature (see Fig. S1 of the supplementary material). However, for shear rate values greater than 2 DPD units an unphysical dependence of the viscosity on the shear rate is obtained. A shear rate of 1 DPD unit is therefore the maximum applicable in our simulation setup. It must be highlighted that the set of parameters used to describe the Pluronic L64 chains is valid in equilibrium conditions and non-equilibrium parametrization may differ such that the predictions of some equilibrium properties could result in unrealistic values. The shear rate in DPD units can be different from the physical shear rate at which an analogous situation is reached. When the value of the shear is greater than 1 DPD unit, the system could be exposed to extreme deformation that leads to non-physical results. To set the lower operating limit, velocity profiles across the simulation box at different shear rates were analyzed. When the shear value imposed on the system is around 10 −3 DPD units, the thermal fluctuations due to the DPD thermostat are masking shear effects and the velocity profile is affected by beads, moving according to the temperature of the system (see Fig. S2 of the supplementary material). This effect was tested on both water and water-Pluronic L64 mixtures; therefore, shear rate values smaller than 10 −3 DPD units cannot be explored. In conclusion, by using a conservative approach, we can set the operating range of shear rates between 0.005 and 1 (see the figures reported in the supplementary material).
Viscosity was obtained by averaging the value of the stress tensor P xy every 100 time steps. All the viscosity values are registered after an initial equilibration phase such that initial fluctuations are filtered. The final value was recorded when fluctuations were around ±0.01 µ DPD , by adjusting the simulation time window. In particular, the trend of the viscosity was recorded during the simulation time and the final value was recorded only when fluctuations were in the order of magnitude of 0.01 µ DPD .
The harmonic potential was finely tuned in order to suppress the formation of over-elongated chains, hence the K Harm coefficient was initially set equal to 4.0 (in DPD units) and then modified. 4 One concentration (i.e., 25% w/w of Pluronic L64 in water) was used as a sample, and K Harm was increased until variations in the viscosity were negligible. In order to reduce potential errors due to extreme shear conditions and over-elongation, the same set of parameters used for the harmonic potential was used in the FENE potential simulations, this means that the value of the spring constant, K FENE , was set equal to 50 DPD units, while the equilibrium distance, r e , was set equal to 1.00 DPD units. Rheograms were obtained for different concentrations of Pluronic L64 at different shear rates, recording the value of viscosity every 0.01 DPD shear units. The qualitative variation of the trend of the viscosity was proven to be related to differences in the microstructure.
Simulations were performed on a cluster InfiniBand 4 TFLOPS on 10 cores AMD Bulldozer and 128 GB of RAM.

IV. RESULTS AND DISCUSSION
In this section, we present the results obtained running DPD simulations for mixtures of water and Pluronic L64 at different concentrations. The first part focuses on the comparison between predictions obtained with equilibrium DPD simulations, on the whole spectrum of Pluronic L64 concentrations, with the experimentally observed structures, 48 whereas the second part shows the effect of shear, investigated with non-equilibrium DPD simulations.

A. Equilibrium simulations
Equilibrium DPD simulations have been performed keeping the temperature bounded at around k B T = 1.00 DPD units (corresponding to 298 K and indicated as a red line in Fig. 1). The snapshots of the DPD simulations are reported in Fig. 1, which also includes the experimentally measured phase diagram. As it can be qualitatively appreciated, peculiar microstructures emerge while concentration increases. Above the critical micellar concentration (CMC), nearly spherical micelles are present at low concentrations (below 25% w/w), while elongated structures can be appreciated at slightly higher concentrations. The presence of a structured network is found when the Pluronic L64 concentration is above 40% w/w, resembling the structure of a soft-gel, while clear lamellar structures are obtained between 70% and 80% w/w. At very high concentrations (i.e., above 85% w/w) small clusters of water are trapped into a Pluronic L64 network, corresponding to the so-called reverse micelles. A closer observation of Fig. 1 confirms that these findings are in good agreement with the measured phase diagram. The length of the box ensures that enough structures (i.e., the number of micelles, hexagons, and layers in the lamellar phase) are obtained to avoid simulation artifacts. As previously discussed, different lengths of the box were tested and an optimal compromise between simulation time and reduction of artifacts was found in the 30 × r c box.
In order to quantitatively characterize the different microstructures, the data collected from the cluster analysis were also used. Figure 2 reports, for example, in the concentration range between 5% and 25% w/w of Pluronic 41 L 1 is for the micellar phase, H is for the hexagonal phase, L α is for the lamellar phase, and L 2 is for the reverse micellar phase. The red line indicates the investigated temperature. Right: the selected snapshots of DPD equilibrium at different Pluronic L64 concentrations (from left to right and top to bottom: 5%, 15%, 25%, 50%, 75%, and 95%). Green beads represent PPO, red beads represent PEO, and white beads represent water and until 40%, only PPO beads are shown. In the last snapshot, only PEO and water beads are shown. As it is possible to see by increasing the Pluronic L64 concentration, larger micelles are formed, as also confirmed in Fig. 3, which also shows that the number of micelles formed increases with the Pluronic L64 concentration. Figure 2 also shows that, at the lowest concentrations and for aggregation numbers greater than four, the micelle radius of gyration scales with the aggregation number with an exponent of d = 3, highlighting the presence of nearly spherical micelles. However, for the largest concentration of the polymer that still allows micellar structures to form (e.g., 25% w/w), the clusters characterized by higher aggregation numbers (above 100) change their geometry from spherical (d = 3) to cylindrical (d = 2), clearly indicating the emerging of elongated micelles.
The data collected from the cluster analysis also confirmed that when the concentration of Pluronic L64 is greater than 40% micelles undergo a transition from spherical to elongated structures, 53 resulting eventually in the formation of a bi-continuous phase. An example is shown in Fig. 4 (left) for a Pluronic L64 concentration of 60% w/w. Finally, at even higher concentrations, lamellae can be observed, but they are not properly identified by the clustering algorithm due to small interconnections between them. More quantitative results from the cluster analysis will be discussed in Sec. IV B, together with the non-equilibrium simulations.

B. Non-equilibrium simulations
As already mentioned, non-equilibrium DPD simulations were used to explore qualitative drawbacks of shear FIG. 4. Snapshots of the DPD simulations for a Pluronic L64 concentration of 60% w/w. An interconnected structure can be observed at equilibrium (left), while the hexagonal phase can be appreciated when the system undergoes shear (right) with a shear rate of 0.1 DPD units.
on the observed microstructures focusing on the variation of the apparent viscosity, number, and dimension of microstructure as a function of the applied shear. Different batches of simulations were run using Harmonic and FENE potentials.
In a DPD model, the parameters for the non-bonded forces are normally selected based on thermodynamic properties such as solubility coefficients or compressibility data. 1 However, a specific criterion to select the bonded parameter (especially when shear is applied) has not been established yet. We decided to tune the value of the harmonic constant, monitoring the viscosity values. A fairly weak spring constant was initially selected, keeping the equilibrium distance equal to 1.00 according to the literature, 49 and then gradually increased. The resulting viscosity values at different shear rates were then recorded and plotted (Fig. 5). When the spring constant exceeded 50 DPD units, further increases had no effect on viscosity. Once this limit was found, before starting non-equilibrium simulations, the FENE potential was used to describe the bonded potential using the same set of parameters obtained from the fine tuning of harmonic potential. This was done in order to reduce over-elongation of the chains when shear values were high. A comparison between results obtained with the two potentials at different Pluronic L64 concentrations and different shear rates is summarized in Fig. 6. As it is clear, small deviation in numerical values is present, but the emerging rheological behavior is similar for both cases. The constants used to describe the FENE potential are K FENE , equal to 50 DPD units, and r e equal to 1.00 DPD units. Now that the bonded parameters have been chosen, we can investigate the effect that the shear rate and polymer concentration have on the viscosity and the microstructure of the liquid. FIG. 3. Observed microstructures at, from left to right, 5%, 10%, 15%, and 25% w/w. Green beads represent the PPO part of the Pluronic L64, while water and PEO beads are not shown. Figure 7 shows how the viscosity changes with the shear rate at different Pluronic L64 concentrations, ranging from 0% (pure water) to 85% w/w. As expected, for pure water no changes in the apparent viscosity are observed, whereas with increasing the Pluronic L64 concentration, the mixture develops a shear-thinning behavior, which becomes more pronounced as the polymer concentration increases. The soft-gel structure, obtained at a concentration of 60%, when the shear rate is applied on the system, is destroyed and an ordered hexagonal phase, constituted by long elongated cylinders perfectly aligned to the streaming flow, emerges. Because of this structural change, a qualitative drop in the viscosity of the fluid can be observed.

C. Deviation from the equilibrium configurations
In this section, we quantify shear effects on microstructures by using the cluster analysis. One example is already reported in Fig. 4, where the formation of the hexagonal phase was qualitatively observed at a Pluronic L64 concentration of 60% by weight under a shear rate of 0.1 DPD units. In Fig. 8, the number of polymer clusters, calculated with the clustering algorithm, is plotted against the simulation time, for three different Pluronic L64 concentrations, representative of three different microstructures, with and without shear acting on the simulation box. The shear rate used in this example is an intermediate value of 0.1 DPD units, far from both the extremely high and extremely low shear regions. Results obtained at other shear rate values led to very similar results. The figure reports simulations obtained with an equilibration phase of 3 × 10 5 DPD time steps, but tests performed with longer equilibration phases (i.e., 2 × 10 6 time steps) and the application of shear for longer times (up to 2 × 10 6 time steps) did not provide relevant differences. As already mentioned, the reported results refer to a size of 30 × r c always leading to the formation of numerous clusters. As it can be seen from Fig. 8, when the Pluronic L64 concentration is around 25% w/w, after an initial transitory phase, in which the randomly positioned FIG. 8. Time evolution of the number of detected clusters in the simulation box for three different Pluronic L64 concentrations of 25% w/w (black), 45% w/w (blue), and 75% w/w (green). Red dashed lines represent the interval in which shear was applied on the simulation box. J. Chem. Phys. 149, 184903 (2018) FIG. 9. DPD simulation snapshots for different Pluronic L64 concentrations (from left to right: 25% w/w, 45% w/w, and 75% w/w) at equilibrium (top) and non-equilibrium (bottom). The shear rate is equal to 0.1 DPD units and only the PPO part of the co-polymer is shown. Different colors represent different clusters found by the clustering algorithm.
chains get closer to each other, spherical micelles are formed, and equilibrium is reached. When shear (represented by the region between the red dashed lines) is applied on the box, its streaming effect induces coalescence between micelles that are close to each other, as proven by the slight reduction in the number of detected micelles/clusters. 54 FIG. 10. Cluster mass distribution (CMD) plotted versus the cluster size or aggregation number detected at equilibrium (top plots in blue) and when shear of 0.1 DPD units is applied (bottom plots in red) for (from top to bottom) Pluronic L64 concentrations of 25% w/w, 45% w/w, and 75% w/w. When the Pluronic L64 concentration is equal to 45% w/w and no shear is applied, a soft-gel is formed and the clustering algorithm detects only one or two clusters. When shear is applied, cylindrical micelles appear, arranged in a hexagonal structure, causing a significant increase in the total number of detected cluster, as clearly visible in Fig. 8 (and in Fig. 4). Similar conclusions can also be drawn for the highest Pluronic L64 concentration (i.e., 75% w/w) for which an increase in the observed number of clusters is visible when the shear is applied.
These modifications in the structures due to shear are also highlighted in Fig. 9, where for three different Pluronic L64 concentrations, the observed cluster is reported at equilibrium (top) and under shear (bottom). As seen for the lowest concentration, the application of shear simply induces micelle coalescence. When the concentration is increased up to 45% w/w at equilibrium, one unique cluster is detected and when shear is applied, the network breaks and hexagonal oriented structures are formed. Finally, the same idea applies to the lamellar phase encountered when the Pluronic L64 concentration is increased up to 75% w/w. At this concentration, the clustering algorithm detects one (or few) cluster, as it is difficult to count the number of lamellas at equilibrium due to interconnections. As visible from Fig. 9, shear is able to break the structure increasing the number of lamellas.
These observations are quantified in Fig. 10 that shows the cluster mass distributions (CMDs), expressed as the frequency of detected structures versus the aggregation number, for the three Pluronic L64 concentrations of 25% w/w, 45% w/w, and 75% w/w with and without shear. As it can be seen, application of shear at the lowest concentration slightly changes the CMD, resulting in the formation of larger (spherical) micelles and increasing the cluster size (or aggregation number) from 30-40 to 60-80. Both at 45% and 75% w/w, the presence of one large cluster or a few larger clusters is detected without shear and the application of shear induces the formation of a few smaller structures.

V. CONCLUSIONS
In this study, we simulated the entire phase diagram of a structured fluid, composed of Pluronic L64 and water by using DPD. We were able to recognize the peculiar structures that form when the concentration varies from dilute to dense, namely, spherical and rod-like micelles, hexagonal and lamellar phases, as well as reverse micelles. Results on this part of the work were found in good agreement with experiments. A novel clustering algorithm was used to identify the structures formed, characterize them in terms of radius of gyration and aggregation number and cluster mass distributions.
Eventually non-equilibrium simulations were also performed, in order to predict how structures are affected by shear, both via qualitative and quantitative analyses. We had to tune bonded interactions between beads, belonging to the same chain, because we noticed that weak springs were affecting the overall behavior of macroscopic dynamic properties when shear was applied. Two different types of potentials were tested and tuned in order to limit this effect. We acknowledge that further investigation is needed to understand differences in the numerical values.
Despite the well-known scaling problem that result in unrealistic shear rates in real units (similar to what happens in non-equilibrium AAMD), we proved that non-Newtonian behaviors can be predicted by DPD and associated with variations of the observed microstructures. Evident drops have been highlighted at higher concentrations, where a transition between phases was more evident.
We are currently investigating different species of Pluronic, using the same set of parameters already used, in order to ensure that the parameters are applicable to different Pluronic.

SUPPLEMENTARY MATERIAL
See supplementary material for the viscosity curve of a system containing water-like beads and the velocity profiles across the simulation box by using Lees-Edwards boundary conditions. Also, the clustering algorithm is explained in this section and can be downloaded through GitHub.

ACKNOWLEDGMENTS
Computational resources were provided by HPC@ POLITO, a project of Academic Computing within the Department of Control and Computer Engineering at the Politecnico di Torino (http://www.hpc.polito.it). The authors want to thank Professor Mike Allen for the useful discussion on scaling factors and a clustering algorithm. One of the authors (I.P.) acknowledges MINECO and DURSI for financial support under Project Nos. FIS2015-67837-P and 2014SGR-922, respectively.