Sarah J.
Gray
,
Martin
Walker
,
Rachel
Hendrikse
and
Mark R.
Wilson
*
Department of Chemistry, Durham University, Lower Mountjoy, Stockton Road, Durham, DH1 3LE, UK. E-mail: mark.wilson@durham.ac.uk
First published on 11th April 2023
Dissipative particle dynamics (DPD) provides a powerful coarse-grained simulation technique for the study of a wide range of soft matter systems. Here, we investigate the transferability of DPD models to the prediction of anionic surfactant phase diagrams, taking advantage of fast parameter sweeps to optimise the choice of DPD parameters for these systems. Parameters are developed which provide a good representation of the phase diagrams of SDS (sodium dodecyl sulfate) and three different isomeric forms of LAS (linear alkylbenzene sulfonates) across an extensive concentration range. A high degree of transferability is seen, with parameters readily transferable to other systems, such as AES (alkyl ether sulfates). Excellent agreement is obtained with experimentally measured quantities, such as the lamellar layer spacing. Isosurfaces are produced from the surfactant head group, from which the second moment M of the isosurface normal distribution is calculated for different phase structures. Lyotropic liquid crystalline phases are characterised by a combination of the eigenvalues of M, radial distribution functions, and visual inspections.
Studying the self-assembly process of surfactants using molecular dynamics (MD) simulations has traditionally been difficult, due to the long time and length scales required. As a result, MD simulations are often performed on pre-assembled surfactant systems,2–7 and studies investigating the full phase diagram are limited. These computational challenges led to the development of the mesoscopic modelling method of dissipative particle dynamics (DPD).8–10 DPD is capable of simulating very large systems while retaining chemical detail, and with its size, speed, and simplicity, DPD has seen many applications since its conception nearly twenty-five years ago, across the broad range of soft matter research.11–18
DPD has been extensively applied to surfactant systems, from simplified dumbbell models,19–21 through to more complex models fitted to reproduce the specific behaviours of individual species.22–25 Within these studies, a wide range of DPD parametrisation methods have been developed. Often these are targeted at dilute solutions where (for example) infinite dilution activity coefficients can be employed,23 or parametrisation can proceed via the surface site interaction point (SSIP) model.26,27
In common with many coarse-grained methodologies, DPD is often found to lack parameter transferability across different thermodynamic conditions, particularly with respect to large changes in concentration and/or temperature. Two-body coarse-grained models (and to a lesser extent – fully atomistic models) are not completely transferable across chemical environments, due to the limits of effective two-body potentials, which average out the effects of many-body forces in the chemical environment used for parametrisation. However, the success of theoretical approaches, such as SAFT28–31 and coarse-grained models such as MARTINI 332–34 make it clear that careful parametrisation can lead to a high degree of transferability. This is particularly true for top–down models parametrised across a range of thermodynamic conditions.35
In this work, we present DPD models of two key industrial surfactants: sodium dodecylsulfate (SDS) and linear alkylbenzene sulfonate (LAS) (Fig. 1).
SDS and LAS molecules have been previously studied at a range of levels from all-atom studies of micellisation36–39 and interfacial behaviour,40,41 to coarse grain studies investigating mesophase structure, energetics, and interfacial properties.42–44 While SDS and LAS have been studied previously using DPD and other coarse-grained methods,43,45–49 only a handful of these previous studies46,47,49 have attempted to reproduce the experimental phase diagram across a range of concentrations. In this paper, we present a parametrised DPD model, which shows a high degree of transferability, not seen in previous works, which performs well across a wide concentration range and readily applies to different molecules. We make extensive use of fast parameter sweeps to optimise DPD parameters to reproduce representative phase diagrams for both SDS and three different isomeric forms of LAS using a set of transferable parameters. The results provide an excellent representation of experimental phase diagrams for SDS and for linear and two-branched forms of LAS, and provide an excellent representation of experimentally measured quantities (phase boundaries and lamellar layer spacings). Moreover, we show the same parameter set can be readily extended to other molecules, through simulations of an isomeric mixture of AES (alkyl ether sulfates).
![]() | (1) |
![]() | (2) |
![]() | (3) |
The random force injects small amounts of energy into the system, and the dissipative force describes the dissipation of energy due to friction – together they act as the ‘DPD thermostat’. In order to recover the Gibbs ensemble this pair of forces is coupled by a fluctuation–dissipation theorem.52 This results in the relations: σ2 = 2γkBT, and ωD = [ωR]2. The DPD thermostat parameters used in this work were γ = 5.625 and σ = 3.354 (using reduced units – see below), which are chosen in line with existing literature.10,45,53 We use the most commonly chosen function for ωR,
ωR = (1 − rij), | (4) |
The defined cutoff of rC = 1 defines the length scale for the system. Since we chose to define the same cutoff for all bead species, our level of coarse-graining is dictated by the smallest group in our surfactants. DPD typically uses reduced units, which require scaling to map to a real system. In this work, masses are scaled such that beads have a reduced mass of m = 1 and energies E are defined by setting kBT = 1. The dimensions of various DPD parameters are provided in Table 1, where various parameters can be converted into real units using their base units.
The bead connectivity and molecular shape are maintained by a combination of two-site and three-site harmonic springs which take the form:
![]() | (5) |
![]() | (6) |
LAS models involve a straightforward extension of our description of SDS. We assume that the closely related sulfate and sulfonate head groups can be represented by the same bead type, with an additional bead inserted between the head and tail beads to describe the benzene ring linker group in LAS. Bonds and angles take the same form as SDS. For branched forms of LAS, two bonds are made between the chains and the benzene linker group, together with associate bond angles. We do not include an additional tail-benzene-tail angle constraint at the branch point. Bonds and angle interactions are specified in full within the ESI.†
The choice of fixing the cut-off distance for all beads rC = 1, effectively approximates the volume of all bead types to be equivalent. This approach is common in similar DPD models23,46 and is reasonable given the similarity of the volumes occupied by different beads. We expect that the approximation has little impact on transitions between liquid crystalline phases, for example between hexagonal and lamellar structures; however, we appreciate that the volume occupied by parts of the surfactant is more important in the micelle regime, for describing factors such as micellar shape, which is not considered in this work.
We primarily establish the undetermined nonbonded interaction parameters for SDS through a comprehensive sweep of parameter space. The model was then extended to other systems, as discussed in detail below. To do that we make use of Flory–Huggins parameters χ in order to demonstrate the transferability of our model. Groot and Warren10 suggest assigning aij interaction parameters by relating χ to the molecular fragments that a DPD bead represents. In the case of our fragments, χ data is not experimentally available. Instead, similarly to the method of Maiti and McGrother,59 we related the theoretical χ parameter of our fragments to their Hansen solubility data,60 taking the form
![]() | (7) |
We calculated our isosurfaces using the utility built into the DL_MESO software,61,62 which we will briefly summarise. The simulation box is overlaid with a regular lattice with spacing 0.25rC, on which the isosurface density will be calculated. The bead volume is smeared according to a Gaussian smearing function
![]() | (8) |
![]() | (9) |
The three eigenvalues ε1, ε2 and ε3 of this symmetric tensor M are indicative of the mesophase present. Note that these eigenvectors are normalised such that their sum is one.
We find that the mesophases can be categorised as follows:
• if ε1 ≃ ε2 ≃ ε3 generates an isotropic phase, such as micellar or bicontinuous cubic;
• if the phase has a hexagonal symmetry, i.e. hexagonal columnar (or inverse) phase;
• if ε1 + ε2 ≤ ε3 we have a lamellar phase.
This measure, along with visual inspection, can aid the identification of the mesophases formed, which is a useful tool in iteratively improving a DPD model to reproduce experimental phase behaviour. Phase structures can additionally be characterised using pair distribution functions g(r), which describe density variation as a function of distance r, which are used in this work for calculating lamellar layer spacing. The g(r) function is computed by calculating the particle–particle distances between tail bead pairs, at every time step. For computational purposes, these distances are ‘binned’ to create a histogram of particle–particle separation distances. The g(r) function is calculated as
![]() | (10) |
![]() | (11) |
The layer spacing which is calculated from the distribution function g(r), can be confirmed by calculating d via another method. Suppose angle θ is the polar angle to the director of the surfactant molecules (where the director is easily calculated as the average direction vector of the molecules). Due to the periodic boundary conditions, the layers must satisfy Lcos
θ = κd, where κ is the number of layers in the simulation box.49 The combination of two approaches to calculating d allows us to be confident in our calculation of the lamellar spacing.
Initial test simulations were performed in the high concentration regime (20% water and 80 wt% surfactant), where SDS exhibits a lamellar phase.63 In this regime, we expect better transferability of DPD parameters across the phase diagram in comparison to low concentrations. The sweep of parameter space was performed in two stages. We began by varying individual cross-interaction parameters aij, while maintaining all other cross interactions at aij = 25.
• Increasingly repulsive values of aTW were found to favour the aggregation of surfactant molecules.
• Increasingly repulsive values of aHT also favoured aggregation but to a lesser extent than aTW, noting the greater configurational freedom of water molecules in comparison to the bonded head beads.
• On its own, aHW has only a small effect on surfactant–water demixing in the higher concentration regime.
Next, we looked at the effect of varying two aij interaction parameters simultaneously to build up a picture of how synergetic effects influence aggregate formation. A sweep through aHT and aTW parameter space (Fig. 3(a)) indicated that aggregation was favoured by simultaneously high values of aHT (>35) and aTW (>40). A visual illustration of the change in aggregation behaviour is shown in Fig. 4. Here, the unusual ‘pseudo-micellar’ behaviour seen for aTW = 45 and aHT = 45 describes a fluid with a self-assembled nanostructure. This corresponds to the formation of very small aggregate nanostructures, as might be seen in a structured tertiary fluid,64 and cannot be classified as true micellar behaviour. In contrast, a sweep through parameter space for aHW and aTW (Fig. 3(b)), favour values of aHW < 20 for phase formation.
![]() | ||
Fig. 3 Schematic representation of the qualitative phase behaviour of our surfactant model. (Simulations conducted using 80 wt% surfactant beads and 20 wt% water). |
![]() | ||
Fig. 4 Impact of varying parameters aTW and aHT on the self-assembly behaviour. Beads are coloured according to their type: water (cyan), alkyl tail groups (purple), and head groups (orange). Cases are categorised to correspond with Fig. 3 as follows: homogeneous (aTW = 30, aHT = 30), phase separation on the bead length scale (aTW = 30, aHT = 45), phase separation (aTW = 45, aHT = 20), pseudo-micellar (aTW = 45, aHT = 45). (Simulations conducted using 80 wt% surfactant beads and 20 wt% water.) |
Having narrowed down the range of phase stability through Fig. 3, the final choice of parameters was made based on a successful prediction of approximate experimental phase boundaries for the hexagonal and lamellar phases. The final set of parameters chosen is shown in Table 2. The lamellar phase that results from a simulation of 80 wt% surfactant is shown in Fig. 5.
a ij | Sulfonate (H) | Alkyl (T) | Water (W) |
---|---|---|---|
Sulfonate (H) | 25 | 50 | 15 |
Alkyl (T) | — | 25 | 45 |
Water (W) | — | — | 25 |
![]() | ||
Fig. 6 (a) A comparison of experimental phase boundaries (from Huddson et al.65 and Kékicheff et al.66) and those assigned via simulation for sodium dodecyl sulfate. Experimentally phases are identified as micellar (L1), hexagonal (H1), a region of mixed phases (62–69%), lamellar (Lα) and crystalline phases (88%+). There is also a narrow two-phase (L1 + H1) region (36–39%). (b) Also plotted is the variation of eigenvalues εi, which can be used to identify the mesophase present. |
The hexagonal phase is found to be more challenging to form than either the micellar or lamellar phases. This is attributed to the fact that the system needs to form large columnar aggregates, and then align these aggregates throughout the entire domain. As these columns can connect across periodic boundary conditions, most often not aligning with perfect hexagonal symmetry, large amounts of energy are required in order to break and reform the columns in the most favourable geometry. As such, a hexagonal columnar phase, that forms spontaneously, requires a larger system and longer simulation time. Therefore, in order to speed up the equilibration process, we chose to apply a gentle shear to the simulation box, which encourages hexagonal phase formation. Shear is then removed and the system is allowed to re-equilibrate in order to produce the equilibrium structure. In addition to the hexagonal phases, we occasionally require the application of shear in order to encourage the formation of lamellar bilayers. For lamellar phases which lie close to the hexagonal-lamellar phase boundary, we sometimes fail to see the formation of periodic layers under equilibrium conditions. In these cases, no clear phase structures appear, even after considerable running time. Following the application of shear, lamellar structures are able to form with layers which are parallel to the sides of the simulation box.
Fig. 8 shows typical g(r) plots for the lamellar phase of SDS, from which the layer spacing can be extracted. Layer spacings are presented in Table 3 as a function of concentration, calculated using g(r). We estimate our simulation length scale based on matching the water density simulated, to an experimentally known water density at 80 °C.48 This produces a length scale of rC ≈ 6.52 Å. This provides a mechanism for converting layer spacings in DPD units into real units. Experimentally SDS solutions at 75 wt% (T = 65 °C) have a reported layer spacing of 3.4 nm.67 This compares excellently with the value at 75 wt% calculated in this work of 3.3 nm. The lamellar phase shows a decrease in spacing as the surfactant concentration increases, which is also in agreement with experiment.67
![]() | ||
Fig. 8 Pair distribution functions g(r) for lamellar phases of SDS with varying concentration. The difference in peak maxima measures the changing average layer separation as concentration is varied. |
wt% | Layer spacing/rcut | wt% | Layer spacing/rcut |
---|---|---|---|
60 | 5.37* | 80 | 5.03 |
65 | 5.44 | 85 | 5.03 |
70 | 5.29 | 90 | 4.80 |
75 | 5.03 |
In summary, we find that the parameters presented in Table 2 reproduce the SDS phase diagram extremely well. This result alone is not surprising, as the parameters of our model were optimised to reproduce SDS in particular. The success of this model lies in its transferability, not just across a full range of surfactant concentrations but also to different molecules.
The phase behaviour of LAS is much more complex than that of SDS, as LAS refers to a whole family of molecules of different chain lengths, degrees of functionalisation, and positional isomers. In order to make a comparison to experimental data, we simulated three specific isomers of LAS, as detailed phase diagrams were unavailable for characterised isomeric mixtures. We simulated an entirely linear isomer, sodium dodecylbenzene-1-sulfonate (DBS1), and two branched isomers, sodium dodecylbenzene-4-sulfonate (DBS4) and sodium dodecylbenzene-6-sulfonate (DBS6). We chose these particular isomers as they cover the broadest range of phase behaviours exhibited by LAS.68
Flory–Huggins parameters χij for all bead pairs are calculated from Hansen solubility parameters60 using eqn (7) and are presented in Table 4. A comparison of the Flory–Huggins parameters for the existing interactions (χHW, χTW and χHT) with the Flory–Huggins for the benzene interactions (χBH, χBT, and χBW) allows us to determine the values of aij for the benzene interactions as aBH = 45, aBT = 25 and aBW = 40. Further information on how additional bead interactions are calculated can be found within the ESI,† including the Hansen solubility parameters used to obtain Table 4.
Sulfonate (H) | Alkyl (T) | Benzene (B) | |
---|---|---|---|
Water (W) | 4.98 | 9.49 | 8.80 |
Sulfonate (H) | — | 3.48 | 2.67 |
Alkyl (T) | — | — | 0.792 |
The DPD phase diagrams of these three molecules are presented (alongside their respective experimental phase information) in Fig. 9. The linear isomer DBS1 shows a similar phase progression (micellar → hexagonal phase → lamellar phase) to that of SDS, while the branched isomers show markedly different behaviour, including the absence of a hexagonal phase. We also note here that all of our micellar solutions (for every isomer simulated) consist of elongated, worm-like micelles as opposed to spherical micelles. This is in agreement with experimental observations of LAS solutions with concentrations 5 + wt%.69
![]() | ||
Fig. 9 Phases assigned to simulations of positional isomers of sodium para-dodecylbenzene sulfonate. Experimental phase boundaries are taken from Ma et al.68 at T = 80 °C. Note that experimental phase region denoted ‘C’ consists of a mixture of liquid crystalline phases, and that ![]() |
The final difference we find between model and experiment is the formation of non-crystalline structure at higher experimental concentrations (90+%). However, this behaviour only occurs in a very small region of the phase diagram at a very high concentration, which is less significant in applications of these models as these conditions are rarely encountered.
At high concentrations 85 + wt%, the experimental phase diagram is complex with many phases present (including lamellar, inverse cubic and inverse hexagonal structures). Here, we do not expect an exact match between experiment and simulation. Primarily, as discussed previously, mesophase behaviour is not on the same length scale as our simulations, and our simulations would struggle to produce coexisting phase structures. Also, as previously discussed, our DPD model is not expected to target phase behaviour in the very high concentration regime.
wt% | DBS1 | DBS4 | DBS6 |
---|---|---|---|
30 | — | — | 13.1 |
35 | — | — | 10.7 |
40 | — | — | 8.93 |
45 | — | — | 7.03 |
50 | — | 6.44* | 7.80 |
55 | — | 6.31 | 6.31 |
60 | 6.44* | 5.98 | 5.69 |
65 | 6.31 | 5.52 | 5.37 |
70 | 6.31 | 5.29 | 5.09 |
75 | 5.98 | 5.37 | 4.85 |
80 | 5.88 | 5.09 | 4.75 |
85 | 5.69 | 4.80 | — |
90 | 5.61 | — | — |
95 | 5.37* | — | — |
We find that for the same concentration, branched molecules have a smaller layer spacing than the linear case, and the DBS6 case has a smaller spacing than DBS4. This is to be expected, given that the end-to-end distance is made smaller with the branching. For all molecules, we find that the lamellar phase exhibits a reduction in layer spacing as surfactant concentration increases. DBS6 solutions continue to form lamellar phases at relatively low concentrations and the spacing grows rapidly with decreasing concentrations. Therefore we only present lamellar spacing values at concentrations 30 + %, since at lower concentrations the spacing is so large that the layers struggle to fit inside the simulation box size. This leads to warping of the layers, which are not perfectly parallel or flat, and the spacing calculation becomes inaccurate. There are also fewer ‘available’ spacing options provided by eqn (11), as the spacing grows. For accurate spacing calculations at lower concentrations, a much larger box size would be required. A visual representation of the results in Table 5 and the available spacings can be found in the ESI.† When the concentration is higher (and the lamellar layer spacing is smaller), there are a significant number of spacings ‘available’, and therefore box size is not a significant issue for this calculation.
It is of interest that the change in lamellar layer spacing for DBS6 is discontinuous, with an anomalous increase at 50% concentration. In the experimental phase diagram, this roughly corresponds to a region in which there are two co-existing lamellar phases (Lα and ). In our simulation, co-existing phases are unable to form (on account of the relatively small simulation length scale). However, it is possible that this is the cause of our anomalous lamellar layer spacing at this concentration. It is worth noting that this discontinuous behaviour is not a result of constraints due to the finite box size (see ESI†).
Existing experimental literature values for lamellar layer spacing of LAS are from isomeric mixtures in water. At a temperature of T = 80° and concentrations 55 + wt%, the lamellar layer spacings are reported as between 3.1 and 3.6 nm.71 By comparison, our data in Table 5 measures layer spacing as between 3.1 and 4.2 nm, in reasonably good agreement with the available experimental data.
![]() | ||
Fig. 10 Coarse-grained representation of alkyl ether sulfates (AES). The number of ethylene oxide (E) beads and alkyl beads (T) are varied in the simulation. |
No. carbon atoms in alkyl tail | No. ethoxy groups | % wt |
---|---|---|
12 | 0 | 29 |
12 | 1 | 14 |
12 | 2 | 7 |
14 | 0 | 29 |
14 | 1 | 14 |
14 | 2 | 7 |
Experimentally, AES exhibits simple phase behaviour, moving from a micellar solution at low concentrations to a hexagonal mesophase at a concentration c ≃ 25%. At higher concentrations, a transition is seen to a lamellar phase at c ≃ 60%.49 The phase progression of AES, as determined using DPD, is described in Fig. 11, and is found to be in quite good agreement with experimental data, with simulations showing the expected L1–H1–Lα phase sequence.
![]() | ||
Fig. 11 Phases assigned to simulations of AES. Experimental phase boundaries are taken from Hendrikse et al.49 at 25 °C (data unavailable at higher temperatures). |
An interesting point raised in this study is that while our models were parameterised via their interaction parameters, the connectivity and shape of the molecule also influence the phase diagram. With no further parametrisation, simply describing the structural differences in two LAS isomers is sufficient for our model to reproduce two distinct sets of phase behaviour. While a simple understanding of the chemical interactions that drive surfactant mesophase formation is important, an equally important role is played by the shape of the molecule.
By simulating mesophases on these length scales, we can gain a further understanding of these complex phases that can be difficult to ascertain experimentally. For instance, the difference in stability of the hexagonal columnar phase between our SDS and our linear LAS models reflects the influence of the surface area to volume ratio in mesophase formation. One obvious limitation of this work is that we do not account for electrostatic interactions; therefore, we are unable to investigate the direct influence of salt concentration on the phase diagram by the simple addition of ions. In our model, the effects of ions are accounted for entirely via the aij parameters. The addition of salt would have an effect on water activity and thus influence the relative aij parameters between solvent and surfactants (and also to some extent the head group self-interactions). Therefore, within our model, for simulations with additional salt, the interaction parameters would need to be recalculated for a new system. Noting, for example, that a reduction in the water–water self-interaction leads to a movement of the L1 – lamellar phase boundary of LAS to lower concentrations, as seen experimentally for LAS on the addition of salt.72
For the consideration of future work, is interesting to note that it may be possible to extend the phase diagrams presented in Fig. 9 with the use of regression algorithms.73 While in this work we have performed individual simulations for each of the LAS monomers at different concentrations, one may be able to infer the phase behaviour of other LAS monomers not studied, without the need to perform additional simulations, through the use of a machine learning technique, such as the quadratic discriminant analysis algorithm (QDA).
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm01641a |
This journal is © The Royal Society of Chemistry 2023 |