Replica exchange dissipative particle dynamics method on threadlike micellar aqueous solutions

Soft matter such as colloids, polymers, gels, and liquid crystals, which have been significantly applied to many industrial fields, has attracted a significant amount of attention as a functional material and has been actively developed since the 1990s. In particular, surfactant molecules are involved in many industrial and engineering applications. This is because the self-assembly of surfactant molecules can spontaneously result in a variety of micelle morphologies, such as spherical micelles, threadlike micelles, and vesicles depending on the concentration of the surfactants, the surfactant molecule structure, the solvent type, and the temperature, among other physical conditions. These polymorphic transitions of surfactant molecules in a solution can drastically influence the properties of functional materials. Therefore, to develop a useful process applicable to the fields of materials chemistry and engineering, it is crucial to predict and control the self-assembly of surfactant molecules.

A molecular simulation is a suitable tool for studying the self-assembly process of surfactant molecules. In fact, many research studies [1–9] on self-assembled structures have been reported through molecular dynamics (MD) simulations. For example, Chu et al [1] investigated the morphologies of giant multi-headed surfactants. Their simulation showed a theoretical guide to controlling the transitions among different morphologies (spherical/cylindrical micelles and vesicles) by tuning the number and topology of the surfactant heads and the surfactant concentration. Wanga [4] also investigated the effects of the spacer group on the self-assembly behaviors of Gemini surfactants. The results of their simulations showed that the self-assembled morphologies transition from spherical micelles and wormlike micelles into vesicles with a decrease in the spacer length.

However, reproducing the most stable structure is difficult for a molecular simulation because systems tend to become trapped in the local minima, particularly within a low-temperature regime. Therefore, to mimic the hydrate crystallization [10–12], a long simulation time, and even a relatively low molecular substance such as water [13–15], is needed. In the case of water, it is an extremely challenging task to gain insight into the solid–liquid phase transitions such as the melting of ice and the freezing of water. The reason for this is that the system should overcome the free energy barrier around the phase transition point, including the heavy calculation costs of the free energies over a wide range of temperatures. Hence, it is considerably more difficult to identify the most stable self-assembled structure in high molecular weight polymers and surfactant solutions because a long simulation time of approximately microseconds and a large number of surfactants and water molecules, namely, 10⁶ water molecules, are needed.

In contrast, to handle greater length and time scales than in a typical MD, a coarse-grained molecular simulation has been widely used. A dissipative particle dynamics (DPD) method [16, 17] is a type of coarse-grained molecular simulation technique, and is a powerful mesoscopic simulation tool enabling the simulation of events occurring within a millisecond timescale and at micrometer length scales. The DPD method has been applied to simulate various self-assembling soft matter systems, including wormlike/threadlike micellar solutions [18, 19], onion-like vesicular solutions [2], and Janus nanoparticle solutions [20–22]. Here, note that the thorough sampling of a system is not assured at a low temperature because of a large complexity with a coarse-grained representation. Nevertheless, the local-minimum problem used in DPD simulations has yet to be reported in detail.

It is known that extended ensemble methods (e.g. a multicanonical ensemble (MUCA) [23, 24], simulated tempering [25, 26], and the replica exchange method (REM) [27, 28]) are effective computational techniques for solving local-minimum problems. Mitsutake et al [29] provides an overview of introduction to the field. In particular, the REM has been used in MD simulations (REMD) [30] to investigate the solid–liquid phase transitions of water in carbon nanotubes [31] and the solid–liquid phase transitions of bulk Lennard-Jones (LJ) fluids [32] because determining the weight factor for the REM is easier than for MUCA. The REM was developed by Hukushima and Nemoto [27, 28] through a Monte Carlo simulation and was applied to MD simulations by Sugita and Okamoto. [30] With this method, systems under different temperature conditions are used simultaneously, and exchanges of such temperatures at certain intervals are applied. Xie et al [33] conducted a Hamiltonian replica exchange Monte Carlo simulation of only a single glucose oxidase molecule, and compared the performance with that of a temperature replica exchange Monte Carlo simulation. The effects of modulating the folding equilibrium of osmolytes on the thermal denaturation of proteins were also investigated using REMD. [34] In addition, a coarse-grained simulation study combined with the REM within only a narrow temperature range has been reported (T  =  0.20–0.60 $k_{\rm B} T$). Ibrahim et al [35] investigated the peptide interactions within a lipid bilayer using a replica exchange on DPD (REDPD) simulations and calculated the potential of the mean force between each peptide within the membrane using the WHAM algorithm. However, as previously described, there have been no studies verifying the practical effectiveness of REDPD when compared to a typical DPD by simulating a long temperature range. In addition, no studies have attempted to apply a larger complex system such as an aqueous surfactant solution, despite the REM being an effective tool for the sampling of the most stable of structures.

In this study, we conducted a REDPD simulation for threadlike micellar aqueous solutions. A threadlike micellar system, which has a wide dynamic range, from the extremely fast motion mode of a single water/surfactant molecule to the thermodynamic fluctuations found in a network structure by multiple threadlike micelles, is extremely complex. We demonstrate that the REM is extremely effective in investigating the most stable self-assembled structures of soft matter when compared to the results of DPD studies without the use of the REM. Moreover, the differences in REMD simulations are discussed.

2.1. Dissipative particle dynamics (DPD) method

We used the dissipative particle dynamics (DPD) method [16, 17], which can simulate millisecond timescales and micrometer length scales to study the self-assembly behavior of a threadlike micellar aqueous solution using an in-house code. The fundamental equation used in the DPD method is Newton's equation of motion. Newton's equation of motion for the particle i is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_i \frac{{\rm d} {{\bf v}_i}}{{\rm d}t} = {{\bf f}_i} = \sum_{j \neq i} {{\bf F}}_{ij}^{\rm C} + \sum_{j \neq i} {{\bf F}}_{ij}^{\rm D} + \sum_{j \neq i} {{\bf F}}_{ij}^{\rm R}\;, \label{eq:eq_motion} \nonumber \end{align} \tag{ 1 }$$

where m is the mass of a particle, ${\bf v}$ is the velocity, ${\bf F}^{{\rm C}}$ is the conservative force, ${\bf F}^{{\rm R}}$ is the pairwise random force, and ${\bf F}^{{\rm D}}$ is the dissipative force. In addition, a_ij is a parameter with information regarding the material used by the DPD method. Moreover, $r_{\mathrm c}$ is the cutoff distance, and ${\bf r}_{ij} = {\bf r}_{j} - {\bf r}_{i}$ and ${\bf n}_{ij} = {\bf r}_{ij} / \left| {\bf r}_{ij} \right|$. Next, the random force (${{\bf F}}_{ij}^{\rm R}$) and the dissipative force (${{\bf F}}_{ij}^{\rm D}$) are defined through the following formulae:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:FR} {{\bf F}}_{ij}^{\rm R} = \left\{\begin{array}{@{}ll@{}} \sigma \omega^{\mathrm R} \left(\left| {\bf r}_{ij}\right| \right) \zeta_{ij} \Delta t^{-1/2} {\bf n}_{ij}, & \left| {\bf r}_{ij} \right| \leqslant r_{\mathrm c} \nonumber \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0,& \left| {\bf r}_{ij} \right| > r_{\mathrm c} \end{array}\right. \nonumber \end{align} \tag{ 2 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:FD} {{\bf F}}_{ij}^{\rm D} = \left\{\begin{array}{@{}ll@{}} - \gamma \omega^{\mathrm D} \left(\left| {\bf r}_{ij} \right| \right) \left({\bf n}_{ij} \cdot {\bf v}_{ij} \right) {\bf n}_{ij}, & \left| {\bf r}_{ij} \right| \leqslant r_{\mathrm c} \nonumber \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0,& \left| {\bf r}_{ij} \right| > r_{\mathrm c} \;\; , \end{array}\right. \nonumber \end{align} \tag{ 3 }$$

where $\sigma$ is the noise parameter, $\gamma$ is the friction parameter, $\zeta_{ij}$ is a random number based on a Gaussian distribution, and ${\bf v}_{ij} = {\bf v}_{j} - {\bf v}_{i}$. The random force (${{\bf F}}_{ij}^{\rm R}$) is equivalent to the frictional force. The temperature is controlled through a combination of the dissipative and random forces. The weighting functions $\omega^{\mathrm R}$ and $\omega^{\mathrm D}$ are related based on the following formula:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:w_func} \omega^{D} \left(r \right) = \left[\omega^{R} \left(r \right) \right]^2 = \left\{\begin{array}{@{}ll@{}} \left[1 - \displaystyle\frac{\left| {\bf r}_{ij} \right|}{r_c}\right]^2, & \left| {\bf r}_{ij} \right| \leqslant r_{\mathrm c} \nonumber \\ \;\;\;\;0,& \left| {\bf r}_{ij} \right| > r_{\mathrm c} \;\; . \end{array}\right. \nonumber \end{align} \tag{ 4 }$$

The noise parameter $\sigma$ and the friction parameter $\gamma$ are linked to each other based on the 'fluctuation–dissipation theorem', which is defined through the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:fd_theory} \sigma^2 = 2 \gamma k_{\rm B} T, \nonumber \end{align} \tag{ 5 }$$

where $k_{\mathrm B}$ is the Boltzmann constant and T is the temperature.

2.2. Replica exchange method (REM)

In the REM simulation, non-interacting independent constant T simulations, called replica, with different T conditions, are applied for the system of concern. The T conditions are exchanged between replicas of neighboring T conditions at certain intervals in accordance with a detailed balance. The probability of an exchange, $P_{{{\rm exc}}}$, between the ith and j th replicas is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{{{\rm exc}}}(i,j) = \min(1, {\rm e}^{-\Delta}), \nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta = (\beta_j - \beta_i)(E_i - E_j), \nonumber \end{align} \tag{ 7 }$$

where $\beta_i$ and $\beta_j$, and E_i and E_j, are the inverse temperatures $1/({k_{\rm B} T}$) and the potential energies of the ith and j th replicas, respectively. Exchanging T conditions enables replicas to randomly walk within a wide space T resulting in an efficient sampling of space E without getting trapped in a local-minimum state.

2.3. Weighted histogram analysis method (WHAM)

The probability distribution of a state with potential energy E in canonical ensemble $P(E;T)$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(E;T) = \frac{1}{Z(T)}n(E){\rm e}^{-\beta E}, \nonumber \end{align} \tag{ 8 }$$

where $n(E)$ is the density of states in space E and $\beta = 1 / (k_{\rm B}T)$. The normalization factor $Z(T)$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z(T) = \int n(E){\rm e}^{-\beta E} {\rm d}E. \nonumber \end{align} \tag{ 9 }$$

Any mechanical or thermodynamic property X is obtained as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left< X \right>_T = \frac{1}{Z(T)} \int X(E) n(E){\rm e}^{-\beta E} {\rm d}E, \nonumber \end{align} \tag{ 10 }$$

where $X(E)$ is an average value of X under a state of E. To compute $n(E)$, the weighted histogram analysis method (WHAM) [31, 36, 37] is used. Specifically, from the REM simulation, accurate histograms $H(E;T)$ in space E under the given conditions T can be recorded, where E corresponds to one of the states that the system can assume and is thus counted once in the histogram $H(E;T)$. Using the obtained histograms $H(E;T)$ and their sum in space E, namely, $h(T) = \sum\nolimits_E H(E;T)$, under the given T conditions, $n(E)$ can be computed by iterating the following two equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle n(E) = \frac{\sum\nolimits_T H(E;T)}{\sum\nolimits_T h(T){\rm e}^{-\beta E + A(T)}}, \nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm e}^{-A(T)} = \sum_E n(E){\rm e}^{-\beta E}, \nonumber \end{align} \tag{ 12 }$$

where $A(T)$ is the Helmholtz-free energy of a micelle system.

2.4. Simulation model and conditions

In this study, we used a surfactant consisting of three beads: one hydrophilic head group (labeled with the letter 'h') and two other hydrophobic tail groups (labeled with the letter 't'). A water (or solvent) particle consists of a single particle (labeled with the letter 'w'). All particles in the surfactant molecule are connected by harmonic springs. The spring force (${\bf F}^{{\rm S}}$) is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:FS} {{\bf F}}_{ij}^{\rm S} = -k_{s} (|{\bf r}_{ij}|-r_{s}){\bf n}_{ij}, \nonumber \end{align} \tag{ 13 }$$

where k_s is the spring constant, and r_s is the equilibrium bond distance. In this study, we adopted ${k}_{s} = 100 k_{{\rm B}}T/r_{\rm c}^2$ and ${r}_{s} = 0.86 {r}_{\mathrm c}$. The length of the surfactant molecule calculated from the bond length distribution and the radial distribution function is approximately 2.5 in a DPD dimensionless unit. The interactions between any two particles in a solution can be described by $a_{{\rm ww}} = a_{{\rm tt}} = a_{{\rm wh}} = 25 k_{{\rm B}}T$, $a_{{\rm ht}} = a_{{\rm wt}} = 70\ k_{{\rm B}}T$, and $a_{{\rm hh}} = 40 k_{{\rm B}}T$. These interaction parameters are related to the solubility parameters. These parameters were an inspiration to the modelling of a short surfactant such as cetyltrimethylammonium bromide (CTAB) as based on a previous study [38]. It is known that this model can realize not only a stable threadlike micelle formation but also similar dynamics (diffusion coefficient of the surfactant molecules) to that for the experiment. These parameters already have been adopted in many studies [39–41]. In addition, their previous studies showed that the micelle is gradually elongated and the threadlike micelle was stabilized by using a moderate repulsive force parameter between hydrophilic head groups ($a_{{\rm hh}} = 40 k_{{\rm B}}T$). The sizes (radius and mass) of single beads are set to the same value regardless of their type [42, 43]. The noise amplitude was set to 3.0, and the friction coefficient was set to 4.5. The volume of the simulation box was 20 $\times$ 20 $\times$ 20 in dimensionless units. The density of the solution ($\rho$) selected was 5.0, and thus the total number of beads in the system was 40 000, of which 4200 particles (c  =  10.5%) were surfactant molecules and the others were water particles. The periodic boundary condition was set in three dimensions.

In our study, we conducted both typical DPD and REDPD simulations of threadlike micellar aqueous solutions to compare the results from these simulation methods. For a typical DPD simulation, we applied one simulation at each temperature, that is, 35 simulations at various values ranging from 0.30 to 7.10 $k_{{\rm B}}T$ (in increments of 0.2) were used. The total run-time was 40 000 $\tau$ in the dimensionless unit (1000 000 step) at each temperature. The initial configuration for the DPD simulation was prepared at random. For the REDPD simulation, we conducted 15 independent REM simulations of a surfactant aqueous solution with 36 replicas. The T ($k_{{\rm B}}T$) conditions were prepared in logarithmic form for each replica, ranging from 0.300 to 0.499, 0.500 to 0.794, 0.800 to 1.094, 1.100 to 1.394, 1.400 to 1.694, 1.700 to 1.996, 2.000 to 2.202, 2.200 to 2.700, 2.700 to 3.191, 3.200 to 3.691, 3.700 to 4.191, 4.200 to 4.698, 4.700 to 5.198, 5.200 to 5.893, and 5.900 to 7.195, respectively. Note that the temperature range was finely set to start at the temperature at which trapping begins under the local-minimum energy states. The procedure used for setting the temperature range is discussed in more detail below. The total run-time was 40 000 $\tau$ (1000 000 step) for each REDPD simulation, and an exchange step of the temperature between nearest-neighbor replicas was conducted every 40 $\tau$. Here, we checked that the systems reached an equilibrium state by tracking the potential energy of the system and the mean aggregation number of self-assembled micelles over time. In addition, the longest relaxation time (${\tau}_c$) for our surfactant model was estimated in the previous study [41]. We set an exchange step of the temperature between nearest-neighbor replicas as 40 $\tau$ (1000 step), to be over the longest relaxation times (${\tau}_c = 20$). We enforced the REM simulations in order of decreasing temperature range, and thus the initial configuration for each REM simulation was the equilibrium structure obtained by the previous temperature range. For a comparison between the results of the DPD and REDPD simulations, we measured the potential energy and mean aggregation number in our system. The potential energy in the system is defined by summing the DPD force between molecule-molecule and molecule-water interactions. Here, we used WHAM to compute the potential energy and mean aggregation number. Usually, reduced units are adopted to report the DPD results. In this simulation, we used reduced units in terms of the cutoff radius $r_{\mathrm c}$, the particle mass m, and the energy $k_{{\rm B}}T$.

We conducted REDPD and DPD simulations of threadlike micellar aqueous solutions. First, the validity of the simulation conditions for the REDPD was assessed. Figure 1 shows the temperature trajectory of the representative replicas. We have made sure of a trip from low to high temperatures and vice versa at least once on each replica simulation (15 independent REM simulations). The histogram of the potential energy distribution was also estimated from a system under an equilibrium state within a low temperature range (0.300 $\leqq$ T $\leqq$ 0.339 $k_{{\rm B}}T$), as shown in figure 2. We confirmed the overlapping of the potential energy distributions between nearest-neighbor replicas. This means that the replicas can move up and down within a certain temperature range. In other words, a temperature exchange is also not applied if an overlapping of the potential energy distributions does not exist. We also calculated the average acceptance ratio for all temperature ranges to ascertain how much the replica exchange was carried out, which was determined to be approximately 56%. Here, we prepared the temperature range in a fine manner to ensure that the acceptance ratio remained at a temperature of below $T = 2.2\, k_{{\rm B}}T$. Thus, it is important to carefully prepare the temperature range near the temperature at which trapping occurs under local-minimum energy states for not only the REMD but also the REDPD.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the potential energy versus temperature for the threadlike micellar aqueous solutions. The pink line shows data obtained using the REDPD simulation with WHAM, and the green circles indicate representative data obtained using a typical DPD simulation. As can be observed, the data of the REDPD match those of the DPD at T $\geqq$ 1.0 $k_{{\rm B}}T$. However, the potential energy for the DPD simulation was higher than that for the REDPD because the system for the DPD simulation becomes trapped in a local-minimum energy state at a lower temperature. Moreover, graphs of the potential energy versus the reduced simulation time, and representative snapshots of the self-assembled structure of the threadlike micellar aqueous solutions, are given in figure 4 to compare the results of the REDPD and DPD simulations. At a relatively high temperature ($T = 2.5 \,k_{{\rm B}}T$), the potential energy and self-assembled structures (figures 4(a)(1) and (a)(2)) between the REDPD and DPD simulations remain nearly unchanged. However, a difference in the time-average potential energy was observed for the REDPD and DPD simulations (see figure 4(b)). We also found that the self-assembled structures clearly exhibited different results in the REDPD (figure 4(b)(1)) and DPD (figure 4(b)(2)) simulations. In general, a self-assembled structures such as threadlike micelles can be observed at low temperature. However, at low temperature ($T = 0.300 k_{{\rm B}}T$), the most stable self-assembled structure was not obtained as the system applied for a typical DPD simulation became trapped under a local-minimum energy state, as shown in figure 4(b)(2). In contrast, for the REDPD simulation (figure 4(b)(1)), the original stable threadlike micelles can be formed. These results indicate the necessity of applying the REM, even for the DPD in which the energy barrier is comparatively lower than that of the MD. Figure 5 shows the temperature dependence of the mean aggregation number for different simulation methods. At T $\geqq$ 5.500 $k_{{\rm B}}T$, this temperature region indicates the random dispersion of surfactant molecules. Their randomly dispersed surfactant molecules begin to aggregate into several isolated spherical micelles at approximately T $\approx$ 5.500 $k_{{\rm B}}T$. The temperature region of 3.000 $\leqq$ T $\leqq$ 5.500 $k_{{\rm B}}T$ corresponds to the growth of the threadlike micelles through collisions between spherical micelles. When the temperature reaches T $\leqq$ 3.000 $k_{{\rm B}}T$, threadlike micelles are completely formed in the system. A marked difference was observed in the results of the mean aggregation number at the same temperature, depending on the simulation methods applied. The correct values of a certain mean aggregation number were not obtained in the DPD simulation at a low temperature (0.300 $\leqq$ T $\leqq$ 3.000 $k_{{\rm B}}T$). This fact is also confirmed based on the self-assembled structure, as shown in Figure 4b(2). In particular, the correct values of a certain mean aggregation number were not obtained in the DPD simulation at a low temperature (0.300 $\leqq$ T $\leqq$ 3.000 $k_{{\rm B}}T$). This fact is also confirmed based on the self-assembled structure, as shown in Figure 4b(2). Thus, to obtain the most stable self-assembled structure in a soft matter system, it is necessary to use the REDPD, even for the DPD in which the energy barrier is comparatively lower than that of the MD. In addition, conducting simulations during the cooling process is also an effective way to explore the most stable state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A REDPD simulation for threadlike micellar aqueous solutions was conducted using NVT ensembles. In this article, through a DPD simulation, we showed that the REM is necessary to obtain the most stable self-assembled structure of soft matter. The potential energy and mean aggregation number in our system were measured to compare the results for both simulation methods. It was found that higher potential energy for a typical DPD simulation was estimated than that for the REDPD because the system becomes trapped in local-minimum energy states at low temperatures. The self-assembled structures showed clearly different results in the REDPD and DPD simulations. For the DPD, threadlike micelles were not observed at low temperature, unlike for the REDPD. At T $\leqq$ 2.000 $k_{{\rm B}}T$, the DPD resulted in incorrect values of the mean aggregation number. Our simulation results demonstrated the effect of applying the REM, even for the DPD in which the energy barrier is comparatively lower than that of the MD, allowing the most stable self-assembled structure in soft matter systems to be obtained. We believe that our study will contribute to the exploration of the most stable self-assembled structures in many soft matter systems, which can be used to form a variety of assemblies.

YK was supported by JSPS KAKENHI Grant No. 18J20653. NA was supported by JSPS KAKENHI Grant No. 17K14610.