Phase Diagram Study of Catanionic Surfactants Using Dissipative Particle Dynamics.

Dissipative particle dynamics (DPD) simulations has been performed to study the phase transition of a mixture of cationic and anionic surfactants in an aqueous solution as a function of the total concentration in water and the relative ratio of surfactants. The impact of the relative difference between the tail lengths of the cationic and anionic surfactants on the phase diagram has been simulated by tuning the number of DPD beads in the simulation model. This research also discusses the impact of the frequently used values of the parameters associated with the harmonic bonds among the bonded DPD beads on the obtained self-assemblies. We find remarkable differences in the resultant self-assemblies based on different choices of harmonic bond parameters. The performed simulations show an enhanced spectrum of self-assemblies with augmented tail lengths and disparate harmonic bond parameters. The obtained self-assemblies are quite unique and can potentially be used in the future for various applications. We also compare the simulation results of the vesicle structures obtained by modeling the electrostatic interaction in the simulation among the charged beads by explicitly introducing charges with a long-range interaction with those obtained by tuning the implicit electrostatic interaction without the long-range interaction. The effects of the chain length of the model and the harmonic bond parameters on the internal density of DPD beads and stress profiles within the vesicles are examined closely. These results are a significant contribution to understanding the stability of the phases and tailoring of the desired vesicles.

Introduction

Surfactants are amphiphilic molecules that consist of hydrophilic and hydrophobic segments. Surfactant is a portmanteau of surface-active agents due to the high activity on the surfaces. They are frequently used in domestic and industrial applications such as cleaning, hygiene, cosmetics, personal care products, fibers, textiles, paints, plastics, pharmaceuticals, food products, petroleum, and so on. These amphiphilic molecules form various aggregates in aqueous solutions, such as micelles, vesicles, bilayer, and other complex self-assemblies. To further enhance the application of these surfactants, it is necessary to study the factors that modify their phase behavior.[1−5] By changing the various external conditions, the aggregates formed by surfactants can be transformed into one another, such as transforming a micelle into a vesicle.[4] The perpetual dependence on these surfactants makes it much more important to study them thoroughly. There is an abundance of papers on surfactants like sodium dodecyl sulfate (SDS).[6−15] These surfactants have been explored by researchers experimentally as well as computationally. This type of aqueous solution, which is obtained by mixing surfactants of different natures and opposite head charges, is known as a solution of a catanionic surfactant. This terminology was given by Jokela[16] in 1985. These catanionic surfactants are also called “ion-pair amphiphiles”. A close inspection of the studies related to ionic surfactants revealed that mixtures of surfactants with opposite head charges could be of more potential use than a solution that only contains one surfactant in water. In 1985, the term used by Jokela et al. for catanionic surfactant was originally specified for salt-free mixtures only. This division of classes shows a very diverse behavior in the natures of the obtained self-assemblies. This distinctive behavior has been the basis of extensive studies done by various researchers.[17−19]

Catanionic surfactants exhibit interesting properties due to the strong electrostatic attraction between the oppositely charged head groups. This surplus electrostatic potential gives rise to extra flexibility, which further induces the formation of exceptional microstructures that could not be produced in the solution of a single ionic surfactant.[20,21] Experimentally, it has been confirmed that these mixtures of surfactants with opposite charges have astoundingly low critical aggregation concentrations compared to that of the single surfactant.[20,21] The surface activeness of these surfactants is supposed to be higher than that of the single surfactant.[22,23] Compared to double-chained surfactants, catanionic self-assemblies are very sensitive to temperature, organic additives, and salts.[20] Hence, certain experiments have been performed to understand the micelle-to-vesicle transition in catanionic mixtures under external conditions.[24,25] The sensitivity of catanionic solutions to the external parameters could also be utilized to develop templates for the synthesis of hollow particles and model membranes.[26,27] The micelle-to-vesicle transition could be utilized in microreactors and drug delivery.[28−30]

A few researchers have worked on a theoretical continuum-based model to elucidate the various vesicle forms obtained.[31,32] The most commonly used computational technique is molecular dynamics (MD). However, MD has its limitation when it comes to complex fluids due to the associated time and length scales, mainly because of the model interaction potential used between the atoms or molecules. Despite all these constraints, many researchers have performed MD simulations of SDS. A research paper by Shelley et al. based on the MD simulation of SDS was published in 1990.[7] The simulation was run on a 182 ps simulation of a 42 monomer SDS micelle. In 1995, a 120 ps simulation of a 60 monomer SDS micelle was performed.[8] After that, many MD simulation works related to surfactants have been published.[9,10,13,33−35] De Vries et al.[36] reported the MD simulation of vesicle formation in water. However, it was computationally difficult to employ this simulation for large systems due to the limitations of the length and time scales. Similarly, for catanionic mixtures, a molecular dynamics simulation would not be a preferable choice due to the restricted time scale and system size. Considering the above shortcomings of MD, coarse-grain MD could be a boon for the researchers interested in the field of soft matter to study the response or behavior of surfactant self-assemblies or catanionic mixtures. In these simulations, the critical micelle concentration (CMC), density distributions within micelles, and the surface adsorption are calculated using either the MARTINI force field or dissipative particle dynamics.[37−42] The results closely match with both the MD simulation results and the experimental values. In 1993, a solution of catanionic surfactants (SDS and dodecyl trimethylammonium bromide (DTAB) in an aqueous solution) in water was studied experimentally by Herrington et al.[43] A ternary (SDS/DTAB/water) phase diagram was obtained as a result of this research. However, experimental characterization could not be used to analyze the kinetics of self-assemblies obtained. Thus, it acts as a barrier in the fine assessment of the complex structures obtained so far. Comparatively, theoretical and computational tools have proven to be radical techniques when it comes to microscopic insights.

Dissipative particle dynamics (DPD) is a competent coarse-grained simulation technique that has been broadly used in the field of soft matter physics.[44−52] In DPD, the clusters of atoms within a molecule are replaced by appropriate single units, called beads. These beads are subjected to conservative, dissipative, and random forces as per the model. The time evolution is calculated using Newton’s equation of motion. The preference for DPD over other coarse-grained techniques relies on the simplicity and reproducibility of different complex phases in soft matter. The DPD technique was first proposed by Hoogerbrugge and Koelmann in 1992[44] and was first implemented by Groot and Madden in 1998[48] to study a block copolymer. Since then, DPD has been applied to a wide variety of soft matter problems.[53−63] A substantial amount of computational work for self-assemblies has been done in general. However, there are comparatively fewer computational studies specifically based on catanionic mixtures.[58]

In the current research, we explored a catanionic system consisting of cationic and anionic surfactants. We simulated the phase diagrams of self-assemblies of the catanionic mixture as a function of the total concentration in water and the ratio between the cationic and anionic surfactants. We also calculated the effect of tail lengths and harmonic bond parameters on the density and stress profile within a vesicle.

This paper is organized as follows. In section , a prologue of the catanionic mixtures is given. Section explains the simulation details, while section provides details of the simulation models. In section , a delineation of the simulation results is given. Section contains various subsections explaining the details of the simulation results and discussions of the differences and similarities of the simulation results for the various models. Section deals with the conclusion and the future perspective.

Simulation Details

DPD is a coarse-grained simulation technique specifically designed for simulating the hydrodynamic behavior of a polymeric system. DPD is a particle-based technique in which atoms and molecules are lumped (coarse-grained) together as a single entity called DPD-beads. This coarse-graining into DPD-beads is the important feature of this technique, as it plays a major role in increasing the computational speed.[44,45] DPD is similar to MD, but it is possible to analyze complex systems, such as the self-assembly of surfactants, with this technique, which is not feasible with classical methods such as molecular dynamics due to limitations associated with the time and length scales. DPD gives rise to a reduced degree of freedom, hence making the simulation computationally efficient and cheaper compared to others. Thus, it enables the analysis of complex systems at larger time and length scales.

To carry out a DPD simulation, the approach is simple. The primary step is to start with Newton’s second law, which is also the first step in molecular dynamics. However, unlike conventional molecular dynamics, here the force acting between a pair of DPD beads i and j consists of three components: the conservative force , the dissipative force , and the random force .[47] The total force experienced by the bead i, denoted , can be written as sum of the three forces as mentioned below:The forms of these forces are discussed in the following paragraphs.

In the case of ions modeled by DPD beads, and depending on the coarse-graining scheme, can be divided in two parts. One part, , contains all the short-ranged interactions, and second part, , contains the long-ranged explicit electrostatic interaction. In this case, the conservative force acting between the beads is given by , whereand is given by eq and explained in section .

Typically, to reduce the computational cost, electrostatic interactions are also modeled by tuning the value of a to implicitly bring in the effective nature of charge. The terminology adopted for the models used in this study is M1, M2, M3, M4, and M5 where ‘M’ stands for model and the following value denotes the model number. In this study, we adopted explicit electrostatic interactions in M4 and M5 and implicit electrostatic interactions in the rest of the models. These models are explained in the subsections of section . For the beads, which do not contain ions, only is used in the models.

The weight function ωD(r) can be chosen arbitrarily but should satisfy the following relations:where kB is the Boltzmann constant. A simple choice for the weight function is given bywhere the prefactor a in eq is the maximum repulsive conservative force between the particles i and j, γ is the friction coefficient in the dissipative force, and σ is the noise amplitude in the random force . In this paper, we use reduced units only. The mass is given by m = 1 for all the DPD beads, and the unit of length is rc, i.e., the cutoff value. While studying the system, we kept σ = 3 and γ = 4.5, which correspond to kBT = 1, and the time step Δt was kept at 0.03 τ (τ is the time unit). We used the LAMMPS[64] simulator to integrate Newton’s equations of motion using the Shardlow algorithm.[65,66] For bonding between the beads of the surfactant, we used the harmonic spring force given bywhere K is spring constant and r0 is the relaxed length. In literature, different values of the spring constant K and the relaxed length r0 have been used.[11,12,47,56] However, the most frequent values used are K = 4 and r0 = 0 and K = 100 and r0 = 0.7, as discussed by Goicochea et al.[56] We performed simulations using both sets of K to understand the difference in the outcomes. The values of a are different in the different types of models used in this study, as explained in section . In section , we also discuss the required details for the chosen values of a in the corresponding model. The size of the simulation box used here is 20 × 20 × 20, and the number density is ρ = 3; thus, the box contains 24 000 DPD beads. Periodic boundary conditions and the NVT ensemble were adopted. The effect of the box size on the aggregate structure was investigated. In this work, at least 105 DPD steps were skipped at each phase point before any statistical quantity related to the phase behavior was calculated. In Figure , we show the kinetics of the evolution of the vesicle phase in time for M1 Case IV, starting from random position of the surfactants. It is clear from the figure that the vesicle had already formed at 7.5 × 104 time steps. We ran this simulation for 5 × 106 time steps, and system remained in the vesicle phase. We also checked energy and pressure fluctuations to test that the system reached the equilibrium state.[14] These observations clearly indicate that the computation time is long enough for the system to achieve an equilibrium state.

Figure 3

Here we show the kinetics of the evolution of the vesicle phase over time, starting from the random position of surfactants for M1 Case IV. At 7.5 × 104 time steps, the vesicle has already formed. We ran this simulation for 5 × 106 time steps, and the system remained in the vesicle phase. In all Cases, all the calculations on these vesicles were started after at least 5 × 105 steps were run.

Additionally, it is clear from the calculated values of radii of gyrations, shown in Table , that the length of box side is 10 times larger than the size of a surfactant. Hence, we do not expect any finite size effects in the results of our simulation. Despite this, we further performed a simulation for the M1 Case II with box size of 30 × 30 × 30 and a number density ρ = 3, containing 81 000 DPD beads, for 106 time steps to check the phase stability (see Figure ) for the total concentration Ct = 0.1 and ratio of surfactants f = 1 (see Figure ). The result is depicted in Figure and clearly indicates that the obtained self-assembly is a vesicle (two vesicles of different sizes) only. The simulation results were tested for different initial states.

Table 2

a

model	LH	LT	Lh	Ln	RA	RC	PH11	Ph1	A1	PH2	Ph2	A2	γ/3	p1	p2	Δ
M1 Case I	1.23	0.87	1.24	2.02	0.87	0.62	0.43	0.69	0.59	0.33	0.53	0.45	0.596	2.5	5.8	3.3
M1 Case II	1.22	0.88	1.23	2.36	1.00	0.62	0.55	1.06	0.88	0.50	0.96	0.80	0.868	2.0	5.7	3.7
M1 Case III	0.86	0.74	0.89	1.61	0.69	0.45	0.24	0.43	0.34	0.24	0.44	0.35	0.183	3.2	6.0	2.8
M1 Case IV	0.86	0.74	0.89	1.88	0.78	0.44	0.32	0.67	0.50	0.31	0.65	0.49	0.354	2.6	5.9	3.3
M2 Case I	1.21	0.87	1.23	2.00	0.86	0.61	0.38	0.68	0.56	0.36	0.64	0.53	0.81	2.5	5.7	3.2
M2 Case II	1.21	0.88	1.22	2.34	0.96	0.61	0.52	0.99	0.84	0.46	0.88	0.74	0.902	1.9	5.7	3.8
M2 Case III	0.86	0.74	0.89	1.61	0.68	0.45	0.25	0.45	0.35	0.22	0.40	0.32	0.414	3.1	5.7	2.6
M2 Case IV	0.86	0.74	0.89	1.89	0.77	0.45	0.32	0.67	0.52	0.28	0.60	0.46	0.357	2.5	5.6	3.1
M3 Case I	1.13	0.86	1.04	1.88	0.83	0.52	0.38	0.68	0.56	0.36	0.64	0.53	0.640	2.5	5.6	3.1
M3 Case II	1.22	0.87	1.23	2.30	0.96	0.62	0.46	0.86	0.71	0.46	0.86	0.71	0.993	2.1	5.7	3.6
M3 Case III	0.81	0.74	0.78	1.56	0.67	0.39	0.23	0.46	0.34	0.22	0.44	0.33	0.731	3.2	5.8	2.6
M3 Case IV	0.86	0.74	0.89	1.86	0.76	0.44	0.28	0.58	0.45	0.27	0.56	0.43	1.004	2.8	5.7	2.9
M4 Case I	1.22	0.87	1.23	2.00	0.87	0.61	0.45	0.73	0.62	0.41	0.66	0.57	0.702	2.4	5.7	3.3
M4 Case II	1.22	0.88	1.22	2.35	0.99	0.61	0.56	1.09	0.91	0.51	0.99	0.83	1.120	2.0	5.8	3.8
M4 Case III	0.86	0.74	0.89	1.61	0.69	0.45	0.25	0.45	0.36	0.24	0.44	0.36	0.164	3.0	5.8	2.8
M4 Case IV	0.86	0.74	0.89	1.87	0.77	0.45	0.30	0.63	0.49	0.28	0.59	0.46	0.646	2.7	5.7	3.0
M5 Case I	1.20	0.86	1.22	1.98	0.87	0.61	0.43	0.70	0.60	0.42	0.69	0.59	0.823	2.5	5.9	3.4
M5 Case II	1.18	0.86	1.21	2.18	0.91	0.60	0.47	0.86	0.75	0.43	0.78	0.68	0.924	1.9	5.7	3.8
M5 Case III	0.86	0.74	0.89	1.61	0.69	0.45	0.26	0.47	0.37	0.24	0.44	0.35	0.663	3.2	5.9	2.7
M5 Case IV	0.86	0.74	0.89	1.87	0.76	0.45	0.31	0.66	0.51	0.28	0.59	0.45	0.929	2.7	5.8	3.1
M3 Case A	1.25	0.88	1.25	2.10	0.92	0.88	0.49	0.82	0.74	0.44	0.74	0.66	0.279	2.3	5.7	3.4
M3 Case B	0.87	0.74	0.87	1.62	0.70	0.68	0.25	0.47	0.38	0.23	0.43	0.34	1.302	2.9	5.8	2.9
M1 Case II(s)	1.22	0.88	1.23	2.36	0.99	0.62	0.52	0.99	0.82	0.50	0.95	0.79	0.740	1.9	5.7	3.8
M2 Case II(s)	1.20	0.87	1.22	2.28	0.97	0.61	0.50	0.94	0.79	0.48	0.90	0.76	1.140	2.5	5.9	3.4
M3 Case II(s)	1.22	0.87	1.23	2.26	0.95	0.61	0.49	0.90	0.75	0.46	0.84	0.70	1.276	2.7	5.9	3.2
M1 Case VI	0.86	0.74	0.89	2.12	0.86	0.44	0.81	1.96	1.51	0.41	0.99	0.76	–0.271	1.9	5.6	3.7
M2 Case V	1.21	0.88	1.21	2.70	1.08	0.60	2.79	6.22	5.34	0.64	1.43	1.23	1.186	1.3	5.6	4.3
M2 Case VI	0.86	0.74	0.89	2.17	0.87	0.45	0.89	2.16	1.68	0.43	1.04	0.80	–0.094	1.9	5.6	3.7
M5 Case V	1.19	0.86	1.21	2.61	1.07	0.60	1.95	4.22	3.59	0.62	1.35	1.15	1.299	1.6	5.7	4.1
M5 Case VI	0.86	0.74	0.89	2.17	0.87	0.45	0.92	2.26	1.79	0.42	1.03	0.82	–0.049	1.9	5.6	3.7

Here, the bond length averages (BLAs) for the HT bond, the TT bond, and the ht bond and the length of H1Tn are shown in the second (LH), third (LT), fourth (Lh), and the fifth (L) columns, respectively. The names of the models, described in Figure , are shown in the first column. In the sixth (RH) and seventh (Rh) columns, the radii of gyration (RG) are shown for the anionic (H1Tn) and catanionic (ht) surfactants, respectively. According to the scheme described in the text, the calculated packing fractions at the two peaks of the head densities in the vesicle phase are also shown. In the eighth column, PH1 stands for the packing fraction for anionic surfactant at the first density peak from the center of the vesicle, as shown in Figures –24. Similarly, in the ninth column, Ph1, refers for the packing fraction for the cationic surfactant at peak one. In the eleventh and twelfth columns, PH2 and Ph2, respectively, refer to the packing fractions at the second peak. In the tenth and thirteenth columns, A1 and A2, respectively, refer for the average packing fractions of the mixture, as explained in the text. In the next column, the surface tension (γ) from eq is shown for the models. In the fifteenth, sixteenth, and seventeenth columns, the position of first peak of heads p1, the position of second peak of heads p2, and the distance between the two peaks Δ, respectively, are shown.

Figure 1

Details of coarse-grain models adopted for this study. These models are categorized by the number of DPD beads and the values of the spring parameters (K and r0; see eq ). The top line of this figure shows the names of the models. Here M1, M2, M3, M4, and M5 stand for model one, model two, model three, model four, and model five, respectively. The interactions among the beads for these models are shown in Table . For every model except cases A and B of model three, i.e., M3 Case:A/B, the cationic surfactant is modeled by two beads. The hydrophobic part is represented by a green bead (tail), the hydrophilic part is represented by a blue bead (head), the head part is represented by h, and the tail part is represented by t. For the anionic surfactant, the surfactants are represented by four, five, or six beads. Red bead represents the hydrophilic head part, which is represented by H, while the black bead represents the hydrophobic tail part, which is represented by T. The value of the spring parameter connecting the beads is shown below each case. For M4 and M5, the charges on the beads are shown by +ve and −ve written beside the beads. M4 has additional counterion beads with charge.

Figure 2

Self-assemblies obtained from the DPD simulation for model one (M1) shown in Figure . The results of each case are mentioned under the title associated with each name, as explained in section and Figure . The ratio of cationic and anionic surfactants is shown along the horizontal axis and is denoted f(Cation:Anion). The total concentration of surfactants in water is shown along the vertical axis and is denoted as Ct. The regions of various self-assemblies are shown by their respective symbols. Indicative images of the micelle, vesicle , and bilayer phases are shown in Figures –14.

Figure 4

To check the stability of the phase, we performed the simulation for 81 000 particles and a box size of 30 × 30 × 30 for model M1 Case II (see Figure ) for 106 time steps with the total concentration Ct = 0.1 and the ratio of surfactants f = 1. We still found the vesicle phase, but here we received two different-sized vesicles.

Table 1

DPD Interaction Parameters a from Equation Used In Simulations of the Cationic and Anionic Surfactant Mixture for All Models Shown in Figure

M1, M4, and M5	h (cation)	t (cation)	H (anion)	T (anion)	W, W+, or W–
h (cation)	25.0	177.8	25.0	100	25.3
t (cation)	177.8	25.0	100	25.0	151.5
H (anion)	25.0	100	25.0	116.5	25.0
T (anion)	100	25.0	116.5	25.0	75.0
W, W+, or W–	25.3	151.5	25.0	75.0	25.0

Figure 6

Disk-like micelles obtained as a result of simulating M1 Case II (H1T4 and K = 4) (see Figure ) at Ct = 0.05 and f = 0.5.

Figure 14

Cross-sectional view perpendicular to the length of a tube obtained as a result of simulating M1 Case IV (H1T4 and K = 100) (see Figure ) at Ct = 0.20 and f = 0.5.

Catanionic Mixture: The Model

The models used in this study are pictorially described in Figure . All the models, as shown in Figure , describe the mixture of the cationic and anionic surfactants. The names of the models are M1, M2, M3, M4, and M5. The values of the repulsive parameter a for DPD interactions among the beads are shown in Table . These models were further categorized into several cases depending upon the harmonic bond parameters (K, r0) given in eq and number of beads involved in the models of the surfactants. Except the M3 Case A/B, the catanionic surfactant is modeled by only two DPD beads, and the anionic surfactant is modeled by four, five, or six DPD beads depending on the category it falls into. Following the publication by Alasiri et al.,[62,67,68] for M3 Case A/B, both types of surfactants in the catanionic mixture are modeled by four beads. In a cationic surfactant, the hydrophobic part is expressed by a green bead (tail), the hydrophilic part is expressed by a blue bead (head), the head part is denoted by h, and the tail part is denoted by t. For anionic surfactants, the red bead indicates the hydrophilic head part, which is represented by H, while the black bead indicates the hydrophobic tail part, which is represented by T. As discussed previously, the models are further categorized for different values of the spring parameters (K, r0) connecting the beads, as shown in eq . This classification is shown in terms of the case number below every model in Figure . In M4 and M5, explicit electrostatic interactions among the charges are implemented, and the charge on the bead is shown by writing +ve or −ve beside the beads. M4 has additional counterions that combine with the water beads to create W+(+ve) or W–(−ve), whereas a neutral water molecule is modeled by a single DPD bead denoted W (neutral).

The motivation for these models is based on papers published at various times for cationic and anionic surfactants.[58,62] In the following subsections, we will focus on the models and the possible interaction parameters involved for each of the models in more detail.

M1 Case I

In view of their many biological applications,[59−61] mixtures of cationic and anionic surfactants, i.e., catanionic mixtures, have attracted the attention of many researchers.[58,69] Li et al.[58] performed a DPD simulation for catanionic mixtures consists of SDS and DTAB. The objective of their study was to examine the factors responsible for changing a micelle into a vesicle. They developed a DPD model for a SDS and DTAB solution in water and studied the phase behavior and the kinetics of vesicle formation by controlling both the amount of salt and the temperature. In these mixtures, the packing parameter plays a key role.[5] Our study also involves this model (M1 Case I) by Li et al.,[58] all the results of which were reproduced before being subsumed in this research.

In the model given by Li et al.[58] the cationic (DTAB) surfactant was modeled using two beads. One bead forms a hydrophilic part with the chemical formula −CH2–N+(CH3)3, which is represented by h and the second bead forms a hydrophobic part with chemical formula −C11H23, which is represented by t. Although DTAB and SDS have carbon chains formed by the same number of carbon atoms, one can see a big difference in their critical micelle concentrations (CMCs). For DTAB, the CMC is 15.6 mM, whereas for SDS it is 8 mM. This difference in the CMC is attributed to the larger size of the head part of a DTAB molecule, which offers more steric repulsion. Comparatively, due to large size of the head in DTAB, the relative size of the tail in comparison to its head is smaller in DTAB than in SDS. Therefore, this DPD model is adjusted by controlling the DPD interaction parameter a and modeling SDS using four beads, where one bead represents hydrophilic part and the remaining three beads represent the hydrophobic part. The variation between the CMCs of the surfactants in the catanionic mixture is reflected here by the differences in their tail lengths.[70] This model of SDS is denoted H1T3. Chemically, H represents −OSO3 and T represents −C4H9.[70] The DPD repulsive interaction parameter, chosen by Li et al.,[58] for the same type of beads (a) is 25. This choice was made on the basis of the simulation done by Groot et al. in 1997,[47] which suggested that a = 25 produces the experimental compressibility of water when the number of water molecules in each DPD bead representing water (N) is one. The interaction parameters between different beads were obtained from the work of Chen et al.[55] Chen et al.[55] evaluated the mixing free energy for the two species (SDS and DTAB) via a Monte Carlo simulation and used this value to calculate the Flory–Huggins parameter χ.[71] The value of a is given by a = a + 3.27χ. Although other schemes are also available for the explicit or implicit inclusion of the electrostatic interactions,[11,12,51,72,73] in this model this effect is not considered. We checked for the effect of the implicit or explicit inclusion of electrostatic interactions in M2 or M4, respectively. For a = 25, as shown by Goicochea et al.,[56] the first peak of g(r), the radial distribution function (RDF), is around 0.85. Goicochea et al. tested the effect of a on the RDF and found that the RDF was mostly unaltered even when the a was switched from 25 to 50. The phase diagram presented in the study from Li et al.[58] qualitatively matches the experimental results showed by Herrington et al.[43] Recently, Debashish et al. published the effect of small-angle neutron scattering (SANS) on the aqueous mixture of SDS and DTAB at a total concentration of 50 mM.[74] Except the multilamellar phase, the result of this study also matches the results of the simulations, as shown in Figure . However, this model could not capture precipitation in these mixtures, which is a limitation of DPD methods. Moreover, precipitation is generally found only in nearly equimolar compositions or in samples below their Kraft temperature, and even this precipitation can sometimes be blocked, which is described in the Blahnik et al.[75]

Generally, the harmonic interaction acting among the beads, as given in eq , is (K, r0) = (4, 0). However, in the literature, many different sets of values for the harmonic potential parameters are used to incorporate the bonding potential among the bonded beads. Statistically, (K, r0) = (4, 0) and (K, r0) = (100, 0.7) are frequently used in DPD simulations.[11,58] These parameters are generally chosen to perpetuate the average distance between the bonded beads. Here, we studied the distribution and the average value of the distance between the beads, and details are shown in Figure and Table for all models. Goicochea et al.[56] also calculated the values for the pressure and surface tension on a plane interface as a function of the spring parameter (K, r0) and found that the values obtained were essentially the same when using the spring parameter values (K, r0) = (4, 0) and (K, r0) = (100, 0.7). Our results, as shown in Figure and Table , for the average bond length are also pretty close to these values. In this study, the obtained self-assemblies involve the plane as well as the curved interfaces. Hence, we are interested in the effect of the spring parameter (K, r0) on the details of a vesicle phase for a given concentration Ct and mixture ratio f(Cation:Anion), as shown in Figure . Part of this research work is dedicated to the study of the vesicle phase found at total concentration Ct = 0.1 and mixture ratio f(Cation:Anion) = 1.

Figure 5

Bond length distributions (BLDs) for M1 Case I and M1 Case III. BLDs for M1 Case II and M1 Case IV are exactly same as those for M1 Case I and M1 Case III, respectively. For other models shown in Figure , the distributions for the bonds were found to be similar; hence, their graphs are not shown. The average bond lengths and radii of gyrations for all the models are shown in Table .

Figure 15

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for all four cases, i.e., Cases I–IV, at (f, Ct) = (1, 0.10). The legends H and T stand for the relative density of the head and the tail beads of an anionic surfactant respectively, relative to total density (ρ(r)). Similarly, the legends h and t stand for the relative densities of the cationic surfactant with respect to the total density (ρ(r)). The legend 9W is used to represent the relative density of the water bead with respect to the total density (ρ). In this figure, we plotted ρ/3 instead of ρ so that the curve would approach a unit value. In all four cases, it is clear that ρ approaches the value of global density (ρ = 3) at the center of a vesicle. As we move along the radius of the vesicle, the local value of total density ρ modulates and again smoothly approaches the global density (ρ = 3) upon reaching outside of the vesicle. Similarly, instead of PN – PT, we plotted to shift the value of stress by one.

Figure 24

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M1 Case II, M2 Case II, M3 Case II, and M4 Case II at (f, Ct) = (1, 0.10). The legend details are same as those explained in the captions of Figures and 22. The profiles, which were simulated for Case II for some initial conditions and did not disappear up to the 106 time step, seem to have a stable phase along with the phases shown in Figures , 20, 22, and 18.

M1 Case II and Case V

As shown in Figure , in M1 Case II, five DPD beads (H1T4), one head bead (H1) and four tail beads (T4), are used to represent an anionic surfactant, whereas in M1 Case V six DPD beads (H1T5) are used to represent an anionic surfactant; in both models, the cationic surfactant is modeled by two beads only. the rest of the details of M1 Case II and M1 Case V are same as those of M1 Case I.

For catanionic mixtures, one could assume that the attraction between the heads of the cationic and anionic surfactants is due to their opposite charges causing the surfactants to experience a subtle change in flexibility, which results in various self-assemblies different from the self-assemblies found in the aqueous solution of a single type surfactant. The value of the attraction depends on the condition, namely, whether the salt is present in the catanionic mixture. We studied the effect of the head interactions by comparing the results of the simulations of M2, M4, and M5 with those of M1. However, the majority of changes undergone by these catanionic mixtures could be attributed to the difference in the effective tail lengths of the cationic and anionic surfactants. The variation between the CMCs of the surfactants in the catanionic mixture is reflected in M1 Case I by the difference in their tail lengths.[70] In M1 Case II and M1 Case V, we were interested in studying the impact of tail length and the parameters associated with harmonic bonds between DPD beads on the phase diagram of a catanionic system. Therefore, we add DPD beads representing the tail to the model by Li et al. and modified it in M1 Case II and M1 Case V.

M1 Case III, Case IV, and Case VI

As shown in Figure , in M1 Case III, the anionic surfactant is modeled by H1T3 DPD beads (one head bead, H1, and three tail beads, T3), whereas in M1 Case IV the anionic surfactant is modeled by H1T4 DPD beads (one head bead, H1, and four tail beads, T4); in both cases, the bonded harmonic interaction is (K, r0) = (100, 0.7). For M1 Case VI, there are six DPD beads for the anionic surfactants (H1T5), and bonded harmonic interaction is (K, r0) = (100, 0.7). The rest of the details of the models are same as those of M1 Case I.

M2 Cases I–VI

As we mentioned in section , catanionic surfactants exhibit interesting properties due to the strong electrostatic attraction between the oppositely charged head groups.[20,21] This additional electrostatic interaction among the surfactants can be accounted in several ways. One way is to directly model the charge interaction, as suggested by Groot et al. and Gonzalez et al.[72,76] This type of interaction scheme was used in M4 and M5. On the other hand, one could implicitly take this interaction into account by tuning the repulsive parameter a among the charged beads, as done by, for example, Mai et al. or Goicochea et al.[11,56] Goicochea et al.[56] used the value of a = 35 among same charges and the value of a = 15 among opposite charges. To include the impact of charge in the results, the simulation was performed for many values between 25 and 35 for the repulsive parameter a among the heads (H) of anionic surfactants as well as among the cationic surfactant heads (h) and to compensate the attractive interaction between H and h for the values between 15 and 25, and the same phase behavior was found. Hence, we show the result only for a = 35 and 15 for additional electrostatic repulsion and attraction, respectively, among the charges.

In M2, the repulsive interaction a among the heads bearing the same charge increased to 35 from 25 and the value of a among beads bearing opposite charges decreased from 25 to 15. The modified interactions are shown in Table . As shown in Figure , for the anionic surfactant, the H1T3 model (one head bead, H1, and three tail beads, T3) is used in Case I and Case III. Similarly, in Case II and Case IV, the anionic surfactant is modeled by H1T4, while in Case V and Case VI the anionic surfactant is modeled by H1T5. The harmonic bond with spring parameters (K, r0) = (4, 0) is used in Case I, Case II, and Case V, while (K, r0) = (100, 0.7) is used for Case III, Case IV, and Case VI.

M3 Cases I– IV, Case A, and Case B

For M3, the DPD interactions a shown in Table are based on the work by Alasiri et al.[62,67,68] Alasiri et al. used the conductor-like screening model for real solvents (COSMO-RS), a quantum mechanical theory,[77,78] to determine the interaction parameters between the beads as the input for DPD. The interaction parameter between the like beads a was chosen to be 25, and that for unlike beads a was calculated using the Flory–Huggins parameter χ. At first, Alasiri et al.[67,68] calculated the screening charge densities (SCDs) around the molecule and the volumes of the molecules by COSMO calculations using Dmol3. After that, COSMO-RS was used to determine the infinite dilution activity γ of the solute. They derived eq , which connected the parameter χ, using Flory–Huggins theory, to find γ and v. Here v is the ratio of molecular volumes of the solute v and v.Using the relation χ = 0.286(a – a) at ρ = 3 and a = 25 used by Groot and Warren,[47] they determined the value of a and found the interfacial tension and CMCs of SDS and DTAB, which were very close to the experimental values. We used the temperature-dependent values of a given in the publication by Alasiri et al.[68] at a temperature of 25 °C for SDS and DTAB. The value of the DPD interaction a for M3 is shown in Table .

They divided both SDS and DTAB into four DPD beads. The tail of the surfactant was modeled by three butyl groups (denoted T3 in Figure ), and the head was composed of [SO4–Na+] for SDS (denoted A in Figure ) and N(CH3)3+Br–] for DTAB (denoted C in Figure ). This model, in our scheme, is denoted M3 Case A and M3 Case B. Spring parameters (K, r0) = (4, 0) for Case A and (K, r0) = (100, 0.7) for Case B are used. The details regarding the number of beads for M3 Cases I–IV are same as those explained for M1 and M2 in previous subsections, and the DPD interaction parameter a is mentioned in Table .

M4 Cases I–IV

In M4 and M5, we consider the explicit charge interaction between the charged beads. As shown in Figure , the charge on the beads is shown by writing +ve or −-ve beside the beads. Following Groot,[72] Gonzalez et al.[76] and Mao et al.,[12] used a smeared charge approach to model electrostatic interactions. As done in this approach, to remove the divergence at the separation r = 0 between the charged beads, they considered the charge distribution on DPD particles of the formwhere λ is the decay length of the charge, which was chosen as 0.25rc to follow Mao et al.[12] The interaction potential between two charged distributions separated by a distance r from center to center is given by[72]where Z is the valency of ion i, Γ = e2/(kBTϵ0ϵrc), e is the electron charge, and β = r/λ. From Groot,[72] Γ = 20.08N1/3.

This relation leads to the electrostatic force Fe between two charge distributions around charged beads i and j.We used particle–particle–particle–mesh P3M to implement the long-range interaction on the periodic boundary condition. The range of the direct interaction in real space was chosen to be 3rc. The potential is shown in Figure .

Figure 19

Density of the counterions for for M3. The density of counterions ρI(r) is very small compared to the total density ρ(r) shown in Figure ; hence, to show the details of counterions, the total density ρ and the densities of the anion (H) and the cation(h) are multiplied and displaced as described in the legends of panels a–d. As shown in panels a–c, the counterion density decays very fast and saturates near the second peak of head density of surfactants outside the vesicle. This decay profile is shown in panel e. The profile was further analyzed and fit with aexp(− bx) + c, where x is the distance from the peak position of the counterion density ρI. This fit is a good match. The fitting parameters are shown in Table . In panel f, the potential action between the beads is shown. In panel f, the red line indicates the electrostatic potential between two point charges of the same sign, Ue[λR] shows the smeared charge potential with the decay length λ = 0.25 (black line). U[0.25R](repulsive) and U[0.25R](attractive) are the total potential due to the conservative part of the DPD, where a = 25, and the electrostatic potential determined by the exponentially smeared charge between same charges and opposite charges, respectively.

The details of the DPD interactions a among the beads for M4 and M5 are shown in Table . The rest of the details related to number of beads in each case are the same as those explained for M1 in section . In this model, we consider the counterions of the heads of the surfactants, denoted by W+ and W–, explicitly in Figure and Table . In Table , we show that the DPD repulsive interaction a for W+ and W– is the same as that of neutral water W, but these beads are subjected to an electrostatic interaction in addition to the DPD repulsive interaction.

M5 Cases I–Case IV

As we briefly explained in the section , one type of catanionic mixtures is salt free. Therefore, in this model, we removed the counterions present in M4 to model salt-free catanionic mixtures. The rest of the details related to the electrostatic interaction and the DPD repulsive interactions are same as those in M4 and are shown in Table . The details on the number of DPD beads and the spring parameters (K, r0) in each case are the same as those from M2 and are explained in section .

Results and Discussion

First, we start to present the results of the global phase diagram obtained from the simulation of M1 Cases I–IV. Results of the simulations are summarized in Figure . The horizontal axis in each Case in Figure shows the ratio of the concentrations of cationic and anionic surfactants and is denoted f(Cation:Anion). The vertical axis shows the total concentration of surfactants in water and is denoted Ct. The self-assemblies of surfactants obtained using these simulations are shown as different symbols, which are defined and displayed below each Case in the subparts of Figure . The indicative images of the different self-assemblies are shown in Figures –14.

Bilayer obtained as a result of simulating M1 Case IV (H1T4 and K = 100) (see Figure ) at Ct = 0.20 and f = 2.

The simulation results for Case I are shown in the top left corner of Figure . It is clear from the figure that we mainly got four types of self-assemblies, which include the micelles, the rod-like micelle, the vesicle, and the bilayer. The micelle phase shown in Case I is not strictly spherical and is closer to disk-like in some cases. We plan to do a detailed study on these in the near future. For Case I, our results are in close aggreement with the results reported by Li et al.[58] We show only the indicative figures of self-assemblies for some cases. As discussed in the follow sections, we also produces these types of self-assemblies in other cases; hence, the number of images shown is optimized.

M1 Case II

The simulation results for Case II are shown in the top right corner of Figure . As mentioned in section and shown in Figure , the anionic surfactant chain length in Case II (H1T4) is longer than the anionic surfactant chain length in Case I (H1T3). The types of self-assemblies in Case II are broadly the same as those in Case I, but there are also interesting differences. The shape of the micelles here is more like a disk, as shown in Figure . Compared with Case I, for Case II, the region associated with the vesicle phase is smaller and mainly replaced by the bilayer and disk-like micellar phases. By definition, at a low value of f anionic surfactants dominate the solution, whereas at a high value of f cationic surfactants dominate the solution. Since the anionic surfactant chain length is longer in Case II than in Case I, the self-assembly of Case II at the ratio f = 6 is the least affected by this chain length increase in the anionic surfactant. Therefore, for all the concentrations, the self-assembly in Case II is the same as that we get in Case I at f = 6. In this region, we get a rodlike micelle at low concentrations (Ct) and a bilayer phase at high concentrations (Ct), which is the same scenario as in Case I. The indicative images associated with the rodlike micelle and bilayer phases are shown in Figures and 9 respectively. However, as the value of f decreases, one can notice the differences in the phase diagram of the self-assemblies found in Case II and Case I.

Figure 7

Rod-like micelles obtained as a result simulating M1 Case II (H1T4 and K = 4) (see Figure ) at Ct = 0.05 and f = 6.

Figure 9

Bilayer phase obtained as a result simulating M1 Case II (H1T4 and K = 4) (see Figure ) at Ct = 0.20 and f = 1.

M1 Case III

The simulation results for M1 Case III are shown in the lower-left corner of Figure below the results of Case I. The model used in Case III varies from the model use in Case I in terms of the spring parameters. As mentioned in section III and shown in Figure , in Case III, the value of the spring constant (K) is 100 and the relaxed distance (r0) is 0.7, whereas in Case I the value of the spring constant (K) is 4 and the relaxed distance (r0) is 0. As mentioned previously, these values of the spring parameters are commonly used in the literature. The effect of the spring parameter values on the DPD simulation was examined by the Goicochea et al.[56] They explored the effect of the surfactants on the surface tension and pressure of the oil–water planar interface using a DPD simulation. They reported the results of the simulation in the space of the spring parameters (K, r0). According to their study, the surface tension and pressure have almost same values for the choice of (K, r0) = (4, 0.0) and (K, r0) = (100, 0.7). This study by Goicochea et al.[56] was devoted to a flat interface and a low concentration. We were curious whether this state of affairs prevailed in a more general situation. To our surprise, the results for Case III are remarkably different from the results found in Case I, as shown in Figure . In Case III, some new types of self-assemblies appeared, which are denoted as tube and complex forms in Figure . Indicative images from Case III are shown in the Figures and 11. In Figure , we show the indicative image of the vesicle phase. This vesicle phase appeared at f = 0.5 and Ct = 0.05 for Case III, whereas for Case I it was a micelle phase. In Figure , a tube-like structure is shown, which appeared at f = 0.5 and Ct = 0.20 in Case III. It was not possible to put some of the self-assemblies found in the simulation into simple categories such as vesicle or bilayer. Hence, we give them a combined name of complex self-assemblies. We plan to explore these complex self-assemblies more closely in our future studies.

Figure 10

Vesicle obtained as a result simulating M1 Case III (H1T3 and K = 100) (see Figure ) at Ct = 0.05 and f = 0.5.

Figure 11

Tube is obtained as a result of simulating M1 Case III (H1T3 and K = 100) (see Figure ) at Ct = 0.20 and f = 0.5. This view is a lateral cross section.

M1 Case IV

The simulation results for Case IV (H1T4) are shown in the lower right corner of the Figure below the results for Case II (H1T4). The model used in Case IV differs from the model use in Case II in terms of the spring parameters, as mentioned in section and shown in Figure . In Case IV, the value of the spring constant (K) is 100 and the relaxed distance (r0) is 0.7, which is the same as that in Case III (H1T3), whereas in Cases I (H1T3) and II (H1T4) the value of spring constant (K) is 4 and the relaxed distance (r0) is 0. The types of self-assemblies obtained in Case IV are broadly the same as those found in Case III. Here we found the tube-like phase, the vesicle phase, the bilayer phase, and the complex phase, which is the same as that in Case III. However, the regions of the various phases in (f, Ct) space are different. The region associated with the bilayer phase of Case III is considerably reduced in favor of other phases in Case IV, as shown in Figure . At Ct = 0.05 and f = 2, we found the disk-like micelle. At small f and low Ct, it seems that the vesicle became oblate in shape. These phases are unique to the Case IV model. We plan to explore these phases in the future. The indicative images of the various phases of Case IV are shown in the Figures –14. Comparing the results of Cases III and IV with those of Cases I and III respectively, it is clear that, in contrast to the result of Goicochea et al.,[56] in our case the use of spring parameters (K, r0) = (4, 0.0) instead of (K, r0) = (100, 0.7) has a large impact on the phase diagram of the self-assemblies.

Figure 12

Vesicle is obtained as a result of simulating M1 Case IV (H1T4 and K = 100) (see Figure ) at Ct = 0.15 and f = 0.5.

The general trend in the phase diagrams discussed above at low concentrations and a ratio of 0.5 is micelle or vesicle. As the ratio f changes from 0.5 to 6, the trend again goes back to rod-like micelles. This is in agreement with the recent experimental results for the SANS data of SDS-DTAB complexes at different SDS-DTAB ratios, as shown by Debashish et al.[74] We can see the types of self-assemblies in detail by going through the phase diagrams in Figure . We can also see that at the ratio f = 1 the general trend at the low concentration Ct = 0.1 is micelle. As we increase the total concentration, the large aggregates that form vesicles to bilayers start to appear, which is in close agreement with experimental result shown by Herrington et al.[43] In all the Cases considered for M1 in our study, one can notice that the vesicle phase is formed at total concentration Ct = 0.1 and ratio f = 1, as shown in Figure . We will see how the details of this vesicle phase change as we introduce changes to our models in the latter sections.

As discussed in section , there are various applications of these self-assemblies in almost all industries.[24−27] For example, micelle to vesicle transitions could be utilized in microreactors and drug delivery.[28−30] There could be many important aspects for the study of all the self-assemblies found here, and we should explore them all closely. However, we focus first on the vesicle phase. It is clear from Figure that a vesicle is obtained at Ct = 0.10 and f = 1.0 in all four Cases. Hence, we chose to focus on Ct = 0.10 and f = 1.0 in the phase space for a more detailed investigation, and the results of the simulation are presented in the Figures –24. In the following subsections, we present the results of the study of the density profiles of the various components that form the vesicle and the profiles of the normal and the transverse stress difference at (f, Ct) = (1, 0.10).

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M1 Case VI at (f, Ct) = (1, 0.10). The rest of the details are the same as those mentioned in the caption of Figure .

Figure 18

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M4 Cases I–IV, at (f, Ct) = (1, 0.10). The counterions are expressed using W+ and W–. The rest of the details are the same as those mentioned in the caption of Figure .

Table 3

Values for the Fitting Parameters in Figure e and the Decay of the Counterions in M4a

M4	a	b	c
Case I(W+)	0.0082	2.3959	0.0530
Case I(W−)	0.0153	1.7621	0.0530
Case II(W+)	0.0079	1.8056	0.0457
Case II(W−)	0.0125	1.5657	0.0456
Case III(W+)	0.0055	1.6344	0.0512
Case III(W−)	0.0136	1.3660	0.0512
Case IV(W+)	0.0101	1.7099	0.0456
Case IV(W−)	0.0138	1.2755	0.0453

See Figure ). The decay profile is fit to ρI(x) = aexp(−bx) + c, where ρI(x) is counterion density at a perpendicular distance x from the surface of the vesicle (see Figure e).

Figure 17

Snapshots of the simulation for M1 Case V (left) and M1 Case VI (right) at the time step of 6 × 105. For M1 Case V, two micelles are formed. One micelle is small, and the other is cylindrical. For M1 Case VI (K = 100), the vesicle is formed, but its shape is prolate spheroid and the chains seems to align along the major axis of the spheroid.

Figure 20

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M5 Cases I–IV at (f, Ct) = (1, 0.10). The rest of the details are the same as those mentioned in the caption of Figure .

Figure 21

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M5 Cases I–IV at (f, Ct) = (1, 0.10). The rest of the details are the same as those mentioned in the caption of Figure .

Figure 22

Total density (ρ) and difference in transverse and normal stress (PT – PN) as a function of distance r from the center of the vesicle for M3 Cases I–IV at (f, Ct) = (1, 0.10). The rest of the details are the same as those mentioned in the caption of Figure .

Packing Fraction

To better understand Figures –24, we analyzed the data generated during the simulation, and they are concisely presented in Table . In the following lines, we explain the contents of Table . The self-assemblies are very sensitive to the packing parameter (P),[5] which is defined as the ratio of the volume (v) of the surfactant molecule to the volume of a cylinder with volume a0lc, where a0 is the projected area of the surfactant head at the interface and lc is the length of the surfactant molecule.

The value of P is supposed to be P ≤ 1/3 for spherical micelles, 1/3