Gaussian charge distributions for incorporation of electrostatic

Article Cite This: J. Chem. Theory Comput. 2019, 15, 4197−4207

pubs.acs.org/JCTC

Gaussian Charge Distributions for Incorporation of Electrostatic Interactions in Dissipative Particle Dynamics: Application to SelfAssembly of Surfactants Hossein Eslami,*,†,‡ Marzieh Khani,† and Florian Müller-Plathe‡ †

Department of Chemistry, College of Sciences, Persian Gulf University, Boushehr 75168, Iran Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Technische Universität Darmstadt, Alarich-Weiss-Str. 8, 64287 Darmstadt, Germany

Downloaded via NOTTINGHAM TRENT UNIV on August 13, 2019 at 08:54:50 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: The point charges are distributed over the soft dissipative particle dynamics (DPD) beads using a Gaussian of tunable width. Screening the Gaussian smeared charge distributions, with wider Gaussians of opposite charge, splits the electrostatic interaction into the real- and the reciprocal-space contributions. This method is validated against model test systems in the literature. The method has also been employed to study self-assembly in solutions of sodium dodecyl sulfate (SDS) in water. The critical micelle concentration (CMC) and the equilibrium concentration of free surfactants, in solutions with SDS concentrations varying from CMC to ≈20 times larger than CMC, are in close agreement with experiment. The microscopic structure of the micelles and the distributions of its hydrophobic/hydrophilic groups and counterions at the water interface are in agreement with experiment. The dynamics of monomer exchange between micelles and solution is examined in terms of the intermittent and continuous correlation functions for the fluctuation of micelle size with time. Two discrete relaxation processes, whose relaxation times differ by 2 orders of magnitude are found. Using the natural DPD time unit, defined in terms of thermal velocity, the relaxation times are an order of magnitude shorter than experimental relaxation times for monomer exchange and establishment of equilibrium between surfactants in the solution and micelles through diffusion of surfactants. However, experimentally comparable relaxation times are obtained by defining the DPD time scale such that the calculated diffusion coefficient of water corresponds to its experimental value. Groot3 addressed the incorporation of electrostatic interactions in DPD, by replacing the point charges on the beads with a linear charge distribution, using the particle−particle particle−mesh (P3M) method4 for solving the field equations. Later, González-Melchor et al.5 employed a Slater-type distribution of charges on the beads and proposed expressions for electrostatic energy and force using the Ewald summation6 approach. A quite similar problem was noticed by Coslovich et al.7 in the theoretical description of polyelectrolytes. It is known experimentally that the swollen polymer coils interpenetrate easily, i.e., their centers of mass may coincide. In order to develop a model of interpenetrating polyions, Coslovich et al.7 distributed the total charge of the polyion over the volume of the coil, using a Gaussian whose width is of the order of its radius of gyration. In this model, called ultrasoft restricted primitive model (URPM),7 the finite Coulomb energy of Gaussians at short distances ensures that coils centers-of-mass still may overlap. The screening properties of the URPM model have recently been examined by Warren et al.,8,9

1. INTRODUCTION The dissipative particle dynamics (DPD) simulation method was originally introduced by Hoogerbrugge and Koelman1 as a mesoscopic simulation tool for modeling soft matter systems over the time- and length-scales not achievable in atomistic simulations. The method was later reformulated by Español and Warren,2 to study the hydrodynamic behavior of complex fluids. In this coarse-grained scheme, a number of atoms are lumped together into a DPD bead, whose motion is governed by three types of forces; a conservative force, a dissipative force, and a stochastic force. Although the DPD method is popular for simulating fluids of a variety of molecular shapes and complexities, technical problems arise when electrostatic interactions between soft DPD beads are to be included. In fact the singularity of the Coulomb potential at short distances causes strong ion pairing of charged beads, since in contrast to models with repulsive hard cores, the repulsion of DPD is not sufficient to offset the Coulombic attraction of opposite charges. On the other hand, the long-range Coulomb interactions are essential for correctly determining the structural and dynamical properties of many systems of practical importance, such as polyelectrolytes, surfactants, lipid bilayers, and colloidal suspensions. In 2003, © 2019 American Chemical Society

Received: February 20, 2019 Published: May 30, 2019 4197

DOI: 10.1021/acs.jctc.9b00174 J. Chem. Theory Comput. 2019, 15, 4197−4207

Article

Journal of Chemical Theory and Computation

where kB is the Boltzmann constant and T is the temperature. The weight function w depends linearly on distance between beads.

through analyzing the asymptotic behavior of the pair correlation functions. Regardless of the method of solving electrostatic interactions, one should adopt a proper charge distribution scheme and calculate long-range Coulomb interactions in periodic systems of finite size. The P3M algorithm of Groot3 for solving Poisson’s equation on the mesh is an efficient scheme for speeding up the calculations, but for the sake of accuracy estimation one needs to invoke the standard Ewald summation as a basis for comparing the results with.4,10 On the other hand, in the Ewald summation method with Slater charge distributions,5 the correct long-range expressions for energy/ force, i.e., the interactions between Slater-type charge distributions and the compensating Gaussians in the reciprocal space, have not been given (see section 2). In this article, the point charges are distributed over the beads using a Gaussian of tunable width. A slight modification of the Ewald summation method is sufficient to calculate energies and forces between DPD beads. Validating the method against known test systems in the literature,3 we have employed it in modeling the aggregation behavior of a charged surfactant, sodium dodecyl sulfate (SDS), in water. It is worth noting that while the DPD method has been widely applied to examine the self-assembly and phase behavior of uncharged systems,11−14 limited attention has been paid to explore characteristics of charged systems, where long-range interactions play the dominant role in determining their physical properties. For example, there are only two reports in the literature focusing on the DPD simulation of micellization of ionic surfactants.15,16 This complexity stems largely from dealing with long-range electrostatic interactions between soft DPD beads. An examination of the validity of the present method in prediction of self-assembly in solutions of SDS in water further supports the capability of DPD with Gaussian smeared charge distributions for modeling complex behaviors of surfactants in water.

rij l o o 1− (rij ≤ rC) o o o rC w D(rij) = m o o o o o0 (rij > rC) n

where rC is the cutoff distance. The conservative potential is expressed as a soft short-range pairwise repulsion, i.e., 2 l o rij yz o 1 ijj o z o o o 2 Aij jjj1 − r zzz (rij ≤ rC) uij(rij) = m C{ k o o o o o o 0 (rij > rC) o n

uC(rij) =

and (2)

where FijD, and FijR are the dissipative and random forces, respectively, between particles i and j, vij and rij are relative velocity and position vectors of particles i and j, respectively, riĵ = rij/rij(rij being the distance between the particles), w(rij) is the distance-dependent weight function, γ is the friction coefficient, δ is the noise strength, θij is a random variable with a Gaussian distribution and unit variance, and Δt is the time step. Español and Warren2 showed that in order to generate the canonical ensemble, the dissipative and random weight functions should be related through D

R

2

2

w (rij) = [w (rij)] and δ = 2γkBT

zizje 2 1 1 = zizjlBkBT 4πε0εr rij rij

(6)

where ε0 = 8.854 187 82 × 10−12 C2/J·m is the vacuum permittivity, εr is the relative permittivity of the medium, and lB = e2/(4πε0εrkBT) is the Bjerrum length. For water (εr = 78.3) at room temperature lB = 0.711 nm, which is slightly longer than the rC values conventionally used in DPD simulations (adjusted according to the degree of coarse graining). An effective tool to overcome the formation of undissociable ion pairs between soft charged DPD beads is to replace the point charges located at the center of beads with smeared charge distributions. This subject has been address by several investigators.3,5,7,8 Diffuse Gaussian smeared charge distributions were first used by Chialvo and Cummings17 for the improvement of the short-range polarization behavior of water. They expressed analytical expressions for the interaction energies of Gaussian charge distributions. Later Jeon et al.18 formulated Ewald summation expressions for long-range electrostatic interaction of Gaussian charge distributions. A similar approach was developed by Coslovich et al.7 for the treatment of electrostatic interactions in polyelecrolytes. In connection with the proposed formalisms by Jeon et al.18 and by Coslovich et al.,7 we introduce and discuss the Ewald summation implementation of electrostatic interactions in DPD, where the charge over bead i is distributed in the form of a Gaussian smeared charge distribution (cocenter with the DPD bead), i.e.,

(1)

FijR = δw R (rij)θijΔt −1/2 riĵ

(5)

where Aij is the repulsion amplitude. For charged beads there is an additional conservative potential, the long-range Coulomb interaction, between charged DPD beads. The Coulomb potential between two point charges i and j, with charges zie and zje (e being the elementary charge) is expressed as

2. THEORY In the DPD simulation method, the DPD beads interact through soft repulsive pairwise potentials as well as local random and dissipative forces.1,2 The dissipative and random forces are expressed as follows. FijD = −γw D(rij)(riĵ ·vij)riĵ

(4)

ij α 2 yz ρi (r ) = ziejjj zzz kπ {

3/2 2 2

e −α r

(7)

In eq 7 the inverse length α tunes the width of the Gaussian. The long-range Coulomb contribution to the potential energy of a pair ij with Gaussian charge distributions is7,17−19 uC(rij) kBT

(3) 4198

=

zizjlB rij

i 2 yz αrijzzz erfjjjj k 2 {

(8) DOI: 10.1021/acs.jctc.9b00174 J. Chem. Theory Comput. 2019, 15, 4197−4207

Article

Journal of Chemical Theory and Computation This potential remains finite at short distances (uC/kBT = 0.798zizjlBα at rij = 0) and converges to the Coulomb potential at long distances. According to Coslovich et al.7 it is possible to calculate the interaction energies of Gaussian charge distributions in the reciprocal space. However, screening the Gaussian charge distributions with wider Gaussians of opposite charges, to conduct part of the calculation effort to the direct space is more efficient (see section 4). In this way it is possible to optimize the splitting parameter to match the simulation box size and to calculate part of interactions in real space and another part in the reciprocal space. Therefore, similar to the standard Ewald summation method, to improve the convergence of eq 8, we assume that every Gaussian smeared distribution (with the inverse width α), is surrounded by a diffuser oppositely charged screening Gaussian (see Figure 1).

N

ρ̃(k) =

ny n zy ji n k = 2π jjjj x + + z zzzz j Lx Ly Lz z { k

U C(s) α = lB ̅ kBT π

N

∑ zi2 i=1

(14)

3. SIMULATIONS Two set of simulations were done. In the first, a test system (a 0.6 M solution of NaCl in water) was simulated to validate the method. The same system has been examined by several investigators3,5,8,9 to validate their methods. For the sake of consistency the same set of DPD interaction parameters as reported by Groot3 has been used for water, Na+, and Cl− (see Table 1). The density of the system was 998 kg/m 3 (corresponding to a reduced density of 3). The temperature was fixed at 298 K (T = 1 in reduced units). The integration time step of the velocity−Verlet integration scheme20 was 0.05 ps. The natural DPD unit of time τ = rC(m/kBT)0.5 = 3 ps, which implies that time step in reduced DPD units is 0.05/τ = 0.017. The second set of simulations were devoted to electrostatic interactions in a system of practical interest; micelle formation in solutions of SDS in water. Here the DPD parameters for all types of interactions were taken from a recent report by Anderson et al.21 (see Table 1). Because of the fact that different sizes, rC, for different DPD beads are used in this force field, SI units (instead of the conventional reduced DPD units) were used in this part of work. The temperature was fixed at 298 K. The integration time step was 0.05 ps (in this case the natural DPD unit of time τ = rC(m/kBT)0.5 ≈ 2.1 ps, calculated in terms of rC and m of solvent, reported in Table 1). Simulations were done for 200 ns to generate equilibrium configurations. After equilibration, simulations were done up to another 200 ns for data collection.

The reciprocal part of the Ewald summation, UC(k), includes the interaction energies of Gaussian distributions and the compensating Gaussians in the reciprocal space and is expressed as7,18,19

k≠0

(13)

where n = (nx, ny, nz) is a vector of integers, Lx, Ly, and Lz are the side lengths of the box, and k = |k|. Note that in the limit of very sharp Gaussian charge distributions (large α), eqs 9 and 11 reduce to the corresponding equations for the direct- and reciprocal-space contributions to the potential energy in the standard Ewald summation method. It is worth mentioning that in the method introduced by González-Melchor et al.,5 where Slater-type functions (simple exponentials) are used for the charge distribution on the DPD beads,5 it is the interaction of point charges with the compensating Gaussians which are calculated in the reciprocal space (similar to the standard Ewald summation). However, one must take into account the interactions between the Slater type functions and the compensating Gaussians in the reciprocal space. The self-energy, UC(s), due to the interaction between the Gaussian charge distributions on the beads with themselves, is overcounted in the reciprocal-space contribution to the potential energy and should be subtracted from the sum of UC(d) and UC(k). This self-energy is expressed as7,18,19

where UC(d) is the direct part of Coulomb interaction and 1 α ̅ = −2 (α + αs−2)1/2 (10)

1 −k2 /4α̅ 2 e |ρ ̃(k)|2 k2

(12)

In eqs 11 and 12, the lattice vectors in Fourier space are expressed as

The total charge of the Gaussian screening charge distribution (with the inverse width αs; αs i ÅÇ k 2 ÑÖ { (9)

∑

j

j=1

Figure 1. Representation of Gaussian charge distributions over DPD beads. Screening of Gaussian charge distributions with wider Gaussians of opposite charge (shown in red) splits the Coulomb interactions into the direct- and reciprocal-space contributions. The compensating Gaussians are shown in blue. All Gaussians are cocentered with the DPD beads.

l U C(k) = B kBT 2V

∑ zie−ik·r

(11)

where V is the volume of the box and the Fourier transformed charge distribution, ρ̃(k), reads as 4199


Article

Journal of Chemical Theory and Computation Table 1. Force Field Parameters for Simulation of Water− NaCl and Water−SDS Solutions21 (A) Nonbonded Interactions reduced DPD units

bead type i

j

(H2O)3 (H2O)2 (H2O)2 (H2O)2 (H2O)2 (H2O)2 CH3 CH3 CH3 CH3 (CH2)2 (CH2)2 (CH2)2 CH2OSO3− CH2OSO3− Na(H2O)2+

Aij

SI units Aij (kJ/mol)

(rC)ij

(rC)ij (nm)

Water−NaCl Solutiona (H2O)3 25.0 1.0000 150.0 Solution of Sodium Dodecyl Sulfate in Water (H2O)2 25.0 1.0000 195.33 CH3 45.0 0.9775 351.60 (CH2)2 45.0 1.0370 351.60 CH2OSO3− 17.9 1.1170 139.86 Na(H2O)2+ 25.0 1.0000 195.33 CH3 24.0 0.9550 187.52 (CH2)2 23.0 1.0145 179.71 CH2OSO3− 28.5 1.0945 222.68 Na(H2O)2+ 45.0 0.9775 351.60 (CH2)2 22.0 1.0740 171.90 CH2OSO3− 28.5 1.1540 222.68 Na(H2O)2+ 45.0 1.0370 351.60 CH2OSO3− 13.3 1.2340 103.90 Na(H2O)2+ 17.9 1.1170 139.86 Na(H2O)2+ 25.0 1.0000 195.33

0.6450 0.5650 0.5523 0.5859 0.6311 0.5650 0.5396 0.5732 0.6184 0.5523 0.6068 0.6520 0.5859 0.6972 0.6311 0.5650

Figure 2. Electrostatic potential and force between two unit charges of like sign, distributed over Gaussians (with inverse width α = 3 nm−1), in water (lB = 0.711 nm). For the sake of comparison, potentials and forces between DPD beads with linear and exponential charge distributions, proposed by Groot3 and González-Melchor et al.,5 as well as those pertaining to the Coulomb’s law are also shown. The full and dashed curves represent electrostatic potentials and forces, respectively. The pink dotted curve represents the short-range conservative DPD force for mentioned beads (rC = 0.646 nm ≈ 0.91lB and Aij/(kBTrC2) = 25.0). At short distances the short-range conservative force does not allow ion pair formation.

1

(B) Bond Parameters, ubond = 2 kij(rij − r0)2 bond i−j CH3−(CH2)2 (CH2)2−(CH2) (CH2)2−CH2OSO3− CH3−(CH2)2

reduced DPD units r0

kij

and rC = 0.646 nm as reported in Table 1 for a solution of ions in water, at r ≈ 0.975rC the analytical Coulomb force between the DPD beads is equal to the short-range repulsive force. At r = rC ≈0.91lB, the Coulomb energy between the Gaussian charge distributions is 1.04kBT, which means that the thermal energy is able to separate oppositely charged ions, and hence, no artificial ion pairing occurs. At large distances, r > 1.5lB, the analytical curve calculated using the present method recovers the Coulomb law. The analytical expressions by Groot3 and González-Melchor et al.,5 however, converge to the exact Coulomb law only at distances above 2.25 and 3.25 rC, respectively. 4.1. Bulk Electrolytes. This is a test system examined by previous investigators3,5,8,9 to validate their methods. The system is a 0.6 M NaCl solution in water. We have simulated 3000 particles (2804 water plus 98 cation and 98 anion beads). Three water molecules are combined into one DPD bead. The density is 998 kg/m3 (corresponding to the reduced density of 3). The temperature is fixed at 298 K (T* = 1 in reduced units). Here the conventional DPD energy and size parameters are used to reduce the parameter space. The real- and reciprocal-space contributions to electrostatic force between cations (or anions) are shown in Figure 3. The cutoffs for truncation of real- and reciprocal-space contributions to electrostatic interaction were 2rC and nC = 4 × 4 × 4 = 64, respectively. The inverse widths for the Gaussian charge distributions and the screening Gaussians were 3 nm−1 (α ≈ 1.94/rC) and 2 nm−1 (αs ≈ 1.3/rC), i.e., α = 1.664 nm−1 = 1.075/rC. For the sake of comparison, we have done all the computations in the reciprocal space by choosing α = αs = 3 nm−1. When α = αs the real space contribution to the Coulomb energy is exactly zero and the whole computational effort is done in the reciprocal space. The calculated electrostatic forces as a function of distance for both cases are shown in Figure 3. For the sake of comparison, the exact theoretical curve is also shown.

SI units r0 (nm)

0.29 150.0 0.1640 0.39 150.0 0.2204 0.59 150.0 0.3336 0.29 150.0 0.1640 (C) Angle Parameters

kij (kJ/mol·nm2) 1172.0 1172.0 1172.0 1172.0

For all types of angles (i−j−k), the parameters (θ0 = 180°) and (kijk = 5 in reduced units, corresponding to 12.50 kJ/mol·rad2) in equation 1 uangle = 2 kijk(θijk − θ0)2 are the same.

The parameters Aij and (rC)ij for Na+−Na+, Na+−Cl−, Cl−−Cl−, Na+−(H2O)3, Cl−−(H2O)3, and (H2O)3−(H2O)3 are the same.3 a

All DPD simulations were performed using our simulation package YASP.22 An atomic Verlet neighbor list was used, which was updated every 10 time steps and neighbors were included if they were closer than 1.1 nm.

4. RESULTS AND DISCUSSION In Figure 2 we have shown exact electrostatic interaction energies, calculated using eq 8 for α = 3 nm−1 (α = 2.13/lB) and forces for two Gaussian smeared unit charges of like sign. For the sake of comparison, energies and forces, calculated from expressions given by Groot3 and González-Melchor et al.,5 which are based on linear and Slater type charge distributions, respectively, are shown in the same figure. Because of the generality of the results in Figure 2, kBT and lB are used as energy and size reduction parameters, respectively. To examine their long-range behaviors, we have shown energies and forces of the Coulomb law, as well. Compared to the charge smearing methods by Groot3 and by GonzálezMelchor et al.,5 in the present method the energies/forces are more repulsive at short distances (r < 1.7lB). Adopting the conventional DPD repulsion parameter, Aij/(kBTrC2) = 25.0 4200


Article


Figure 4. Radial distribution functions for cation−cation, cation− anion, and water−water beads in a 0.6 M solution of NaCl in water at T = 298 K. The curves and markers represent our results and calculations by Groot,3 respectively. The DPD interaction parameters are reported in Table 1.

Figure 3. Calculated electrostatic force between cation−cation (anion−anion) with Gaussian smeared charge distributions, in a 0.6 M solution of NaCl in water, as a function of distance. The inverse widths of the smeared and screening Gaussians are shown in the figure’s legend. The black curves indicate calculations using the present method, in which a Gaussian of inverse width α = 3 nm−1 is screened by a wider Gaussian of inverse width αs = 2 nm−1. Here part of calculations (within a cutoff distance 2rC) is done in the real space. The inset shows convergence of the energy in the direct space. The red curve stands for the case in which the whole computational effort is done in the reciprocal space. In both cases a number of 4 × 4 × 4 vectors are used in the reciprocal space.

Figure 2), compared to Groot’s method.3 Identical with the findings of Groot,3 our calculated radial distribution functions satisfy the relation gsolvent−solvent (r ) = (gcation−anion (r )gcation−cation(r ))1/2

(15)

Because of using the same short-range conservative interaction parameters for solvent and ions, here gcation−cation(r) = ganion−anion(r) and gcation/anion−solvent(r) = gsolvent−solvent(r). Having validated the method against the known literature, in the next section we apply it to simulation of self-assembly in solutions of a charged surfactant, SDS, in water. 4.2. Micellization of Sodium Dodecyl Sulfate. In this section the micellization of SDS is studied at a constant temperature (T = 298 K), constant volume, and a constant concentration of SDS in water. The mapping scheme from atoms to DPD beads is depicted in Figure 5. In this mapping, one to three heavy atoms (C and O) are lumped together to construct a DPD bead.21 The only exception is the headgroup bead, CH2OSO3−, which consists of six heavy atoms. For water, a combination of two water molecules constitutes a single DPD bead. The Na+ DPD bead is hydrated; it consists of Na+ plus two water molecules. The DPD parameters, tuned against experimental partition coefficients, solubilities, and densities of small molecules, are reported in a recent publication by Anderson et al.21 In this parametrization, the density of water is set to 998 kg/m3 (ρ = 3 in reduced units). Because of the fact that the critical micelle concentrations (CMCs) obtained using simulation methods normally differ (by factors as large as 10) from experiment,24 a wide range of concentrations of SDS in water were used. In each sample a number of SDS molecules were randomly distributed in a cubic simulation box of size 27.5 nm. In all samples there are a total (water plus SDS) of 346 000 DPD beads in the simulation box. The number of SDS and water molecules was adjusted to achieve a prespecified concentration, ranging from 3 to 220 mM. The most concentrated sample, 220 mM, contains 2754 SDS molecules (22 032 DPD beads) plus 323 968 water beads. Using the error estimate expressions for the standard Ewald summation method23 a total number of 7

Compared to the case where all the computations are done in the reciprocal space, screening the Gaussian charge distributions (of inverse width α) with wider Gaussians (of inverse width αs) to conduct part of the computations in the direct space produces results in a better agreement with the exact theoretical curve. Obviously to increase the agreement between the results generated using the former case and the theoretical curve, a larger cutoff in the reciprocal space is needed. In the inset of Figure 3 we have shown the real-space contribution to the electrostatic energy, which converges to zero at r < 2rC. Compared to the Groot3 and GonzálezMelchor et al.5 methods, here the real-space contribution to the electrostatic energy is shorter range, making the computations in the real space less expensive. It is worth mentioning that using the error estimation expressions for the standard Ewald summation method,23 the deviation in energy in the reciprocal space for α = 1.664 nm−1 = 1.075/rC is of the order of 10−3. Reasonably, using a larger cutoff for the truncation of reciprocal-space contribution to electrostatic interaction, nC, leads to a smaller error. However, as Figure 3 shows, the analytical and calculated curves mostly deviate at short distances (r < rC), where the short-range conservative force of DPD is much stronger than the Coulombic force (see Figure 2). In other words, an increase in nC does not lead to a substantial improvement in the accuracy of computations in the reciprocal space. We have shown in Figure 4 the radial distribution functions, g(r), for cation−cation, cation−anion, and solvent−solvent pairs. The g(r) curves in Figure 4 are also compared with the corresponding distributions generated by Groot3 for the same system. Our calculated g(r) for cation−anion pair is slightly more structured than that of Groot’s. This is due to larger electrostatic interactions at short distances in our method (see 4201


Article


Figure 5. Mapping scheme form all atom to DPD beads,21 adopted in this work.

× 7 × 7 = 343 vectors (α ≈ 3 nm−1 and αs ≈ 1.4 nm−1, corresponding to α = 1.27 nm−1) are needed to adjust the precision for reciprocal-space contribution to the electrostatic energy to that for NaCl test system (see section 4.1). The cutoff distance for truncation of real space contribution to electrostatic energy was 2 nm. Long DPD simulations (200 ns) were performed to achieve equilibration, monitored as the stage at which the concentration of monomers in the solution becomes constant. The aggregates are distinguished from monomers based on a distance criterion, i.e., two monomers are regarded to be involved in the same aggregate if their centers-of-mass locate within a distance 1.35 nm from each other. Similar to previous studies in the literature,15,16 the CMC is defined as the lowest concentration of surfactant in water at which first micelles appear in the system. The probability distribution for micelle size, obtained by averaging the number of aggregates of size n, is plotted in Figure 6. At high SDS concentrations, larger micelles are

trations (as large as 50 times the CMC). It is worth mentioning that the mean aggregation size calculated in simulations are typically smaller than those observed in experiment. The experimental25,26 mean aggregation number for SDS is about 60. Aggregates of this size are seen in our simulation box only at high concentration, 20 times higher than the CMC (see Figure 6). We have shown in Figure 7 the concentration of free (n ≤ 3) surfactants, CF, versus the total concentration of SDS, CT, in

Figure 7. Dependence of free monomer concentration, CF, on the total concentration, CT, of SDS molecules at T = 298 K. The open and filled circles represent our calculations and the experimental data,27 respectively. The curve represents the best fit through calculated points, according to eq 16. The open diamonds and squares represent DPD simulation results of Mao et al.15 and Anderson et al.,16 respectively. The dashed line is the bisector of CF and CT axes. The intersection of the dashed line with the full curve indicates the CMC (10.9 mM). Figure 6. Probability distribution, P(n), for finding SDS micelles of size n at different concentrations, shown in the figure’s legend as multiples of the critical micelle concentration (10.9 mM).

the solution. The CMC is known as the highest total concentration of surfactant at which CF = CT. Below the CMC no micelle forms (CF = CT). The intersection point between the curve connecting CF vs CT values in micelle forming systems and the bisector of CF vs CT corresponds to the CMC. For nonionic surfactants the CMC is characterized by the occurrence of a plateau in the graph of CF vs CT. For ionic surfactants, on the other hand, following the initial rise, the decrease in CF identifies the CMC. This turnover in the concentration of free surfactants, at CT > CMC, is due to the Coulomb interactions between free surfactants and the micelles. The trend observed here for the dependence of CF on CT is in agreement with experimental observations27 and DPD simulation results by Mao et al.15 and by Anderson et al.16

formed. At lower concentrations (above the CMC), on the other hand, small size aggregates and micelles, separated from each other by a region of population depletion in Figure 6, are found. The calculated CMC is 10.9 mM, which is 33% higher than the experimental value (8.2 mM).25 The aggregation number (the average number of surfactants in a micelle), over the concentration range examined in this work, shows a nearly linear concentration dependence. Such a linear relationship between aggregation number and concentration for SDS has been found in previous simulations in the literature.16,24 The aggregation number-concentration curve in experiment,25 however, is a nonlinear (concave) curve; the aggregate size increases with concentration and saturates at high concen4202


Article

Journal of Chemical Theory and Computation Another estimate used in the literature to find the CMC from the calculated CF and CT values is through the empirical relation by Sanders et al.,28 i.e., log(C F) = (1 + β)log(CMC) É ÄÅ ÅÅ (1 − β)(C T − C F) + C F ÑÑÑ ÑÑ Å Å − β logÅÅ ÑÑ ÑÑÖ ÅÅÇ 1 − VC ̅ T (16) where V̅ is the molar volume of surfactant (assuming to correspond to a density of 1000 kg/m3) and β is the degree of association of counterions with the micelle forming surfactants. We left both β and CMC as the fitting parameters to fit our calculated CF vs CT points according to eq 16. We obtained β = 0.75 and CMC = 10.9 mM. Our calculated value of β is very close to the reports in the literature for SDS.15,16,28 We have further shown in Figure 7 the dependence of CF on CT, calculated in previous DPD simulations15,16 and in experiment.25 Compared to the calculation by Mao et al.15 (CMC = 19 mM), our calculated CMC value (10.9 mM) is in a better agreement with experiment25 (8.2 mM). Our calculated CMC is, however, higher than the reported values by Anderson et al.16 (6.7 mM). This can be explained in terms of their definition of free surfactants (monomers), compared to ours (n ≤ 3). Expectedly, there exists a lower concentration of monomers than aggregates of size n ≤ 3 in the system. 4.3. Micelle Structure. We have shown in Figure 8 a snapshot of a configuration of a micelle, consisting 50 SDS

Figure 9. Density profiles for different DPD beads, shown in the figure’s legend, as a function of distance from the center-of-mass of the micelles with aggregation numbers ranging between 55 to 60 SDS molecules. The total concentration of SDS in water is 20 times higher than the CMC.

calculations for micelles of sizes between 55 to 65, i.e., those whose sizes closely match with the experimental micelle size (≈60).25,26 The terminal CH3 beads are mostly concentrated in the interior core of the micelle. However, in agreement with experiment,29 the distribution for the CH3 tail group has a tail extending to the region of the sulfur distribution. The hydrocarbon tail constructs a condensed liquid-like hydrophobic core. The outermost part of the micelle consists of the charged hydrophilic CH2OSO3− beads. A comparison of the density profile peaks for hadrocarbon tail and water shows that water penetrates to a depth of ≈0.4 nm from the surface of the hydrocarbon core of the micelle, i.e., the surface of the hydrocarbon core is wet. This interphase thickness closely correlates with the measured experimental X-ray scattering data26 on SDS and small angle neutron scattering data on lithium dodecyl sulfate.30 Compared to the density profile peak for the CH2OSO3−, the position of maximum in the Na+ density profile locates at a larger distance from the center-ofmass of micelle. The density profiles in Figure 9 are in agreement with experimental reports on the radius of hydrocarbon core (1.7 nm) and the radius of micelle (2.23 nm) measured from X-ray scattering.29 Our findings also support experimental observations revealing that the terminal CH3 groups have a broad distribution throughout the micelle.29 A comparison of the density profiles for water, CH2OSO3−, and C show that water molecules are distributed at the micelle surface; the polar CH2OSO3− groups are hydrated and water slightly penetrates to the hydrocarbon shell. More detailed analysis of the structure of micelle at the water interface is performed by examining the radial distribution function, g(r), for Na+−CH2OSO3− pair (Figure 10). All g(r) curves in Figure 10 are calculated for pairs falling into cones of opening angle 60° (whose vertices coincide with the considered beads in the micelles) orienting along a line connecting the center-of-mass of micelle to the considered bead in the micelle. The sharp g(r) peak at 0.45 nm is an evidence for the presence of contact Na+−CH2OSO3− ion pairs. The broad shoulder between ≈0.8 and ≈2 nm represents the second Na+ shell, separated from the first shell by a layer of water (see the g(r) curve in the inset of Figure 10 for

Figure 8. Snapshot of a SDS micelle configuration. Colors are the same as those in Figure 5. For the sake of clarity, water molecules are not shown.

monomers. The terminal CH3 groups are mostly located in the center of micelle, extending to its surface. The charged CH2OSO3− groups are covering the surface of the hydrocarbon core of the micelle and Na+ ions are mostly located at the interface of micelle with water. For the sake of clarity, water molecules are not shown in Figure 8. More detailed information regarding the distribution of different DPD beads in micelle is shown in Figure 9 in terms of the density profiles of DPD beads as a function of distance from the center-of-mass of micelle. In order to be able to compare the results with experiment, we have done 4203


Article


otherwise. The function H(t), however, is unity when a molecule is continuously (from time t = 0 to t) involved in the same micelle and zero otherwise. The function Cn(t) measures the probability that a molecule involved in a micelle of size n at time t was involved in the same micelle at time t = 0, irrespective of the possible breakings in the interim time. The function Sn(t), on the other hand, measures the probability that a molecule involved in a micelle of size n at time t = 0 continuously remains involved in the same micelle until time t. To investigate the decorrelation behaviors of functions Cn(t) and Sn(t), we have examined the exchange of surfactants between micelles and solution during the consecutively saved configurations. As the decorrelation depends on the time interval between the saved configurations, we wrote a high resolution trajectory (with the time resolution of 3 ps) for this purpose. The intermittent correlation function, plotted in Figure 11, for the exchange of monomers between the micelles

Figure 10. Radial distribution function for the Na+−CH2OSO3− pair. (inset) Radial distribution functions for water and the micelle beads shown in the figure’s legend (see the numbering of beads in Figure 5). All the distribution functions are for solution of 220 mM SDS in water at 298 K. The g(r) curves are calculated for ij pairs falling into cones of opening angle 60° orienting along a line connecting the center-ofmass of micelle to the considered bead. The cone vortex coincides with the considered beads in the micelle.

CH2OSO3− water pair with a maximum at ≈0.6 nm). In the inset of Figure 10 we have shown the radial distribution functions for water−CH2OSO3−, water−terminal CH3, water− (C2/C3) bead, and water−water−(C4/C5) beads. Expectedly, the CH2OSO3− head groups are hydrated. Moreover, in complete agreement with experiment, the surface of the hydrocarbon core (C2/C3) of micelle is wet. The depth of water penetration to the hydrophobic core is not enough to wet the C4/C5. A small fraction of terminal CH3 groups is located at the interface of water. The thickest interphase belongs to the spatial distribution of sodium ions at the surface of micelle (2 nm). 4.4. Monomer Exchange Kinetics. The most accepted mechanism for the variation of the size of micellar aggregates is through stepwise addition/removal of molecules to the micelles. This mechanism is supported by simulations and by recent ultrasonic relaxation experiments27,31 on kinetics of micelle formation, in which it is shown that experimental data for kinetic of micelle formation are well fitted by theoretical models based on surfactant addition/removal to/from micelles. In equilibrium, the micellar aggregate sizes are fluctuating around a constant value. The dynamics of surfactant exchange between a micelle of size n and the surrounding surfactants in solution can be analyzed in terms of the following two correlation functions, originally introduced by Rapaport32 for investigation of the hydrogen bond dynamics in water. The intermittent, Cn(t), and the continuous, Sn(t), correlation functions for micelles of size n, are defined as Cn(t ) =

⟨h(0)h(t )⟩ ⟨h(0)⟩

(17)

⟨h(0)H(t )⟩ ⟨h(0)⟩

(18)

Figure 11. Decay of the intermittent correlation function, Cn(t), for exchange of monomers between micelles of aggregation number n and solution. The sizes of micelles are shown in the figure’s legend.

and solution, shows an initial fast followed by a much slower decay regime. The short-time decay does not depend on the micelle size. Over this short time scale secession of surface molecules from the micelles is the dominant process. The slow decay regime is connected to the diffusion. At longer time scales diffusion brings back some of the dissociated surfactants to the micelle. In agreement with experiment, our results show than only a small fraction (≈2%) of dissociation/association events involve the coalescence/breakup into dimmers, trimers, and higher-order aggregates. The calculated relaxation times for the indicated intermittent correlation functions in Figure 11 are 13, 20, and 25 ns for clusters of size 25, 40, and 55, respectively. This relaxation time (of the order of 10 ns) can be compared with the experimental27 relaxation times for monomer exchange equilibrium, obtained from broadband ultrasonic spectrometry measurements (≈200 ns for a SDS solution of the same concentration). In these experiments, the micellar solutions are subjected to a rapid perturbation and the shift in the aggregation number of micelles from their equilibrium values are measured as a function of time. According to the fluctuation−dissipation theorem, the response of a system to an external perturbation is connected to the internal

and Sn(t ) =

where the binary function h(t) is unity when a particular molecule is involved in a micelle of size n at time t, and zero, 4204


Article

Journal of Chemical Theory and Computation fluctuations existing in an equilibrium system. Here, fluctuations in the micelle size in thermal equilibrium are characterized in terms of the correlation function, Cn(t). The calculated relaxation time for the decorrelation of function Cn(t) is shorter than the experimentally measured27 relaxation time by an order of magnitude. This is a consequence of the well-known overall acceleration of a DPD model over an atomistic model or experiment.2,33 We have also plotted in Figure 12 the decorrelation behavior of the function Sn(t). Compared to Cn(t), the function Sn(t)

To have an estimate of the magnitude of dynamic acceleration in our DPD simulations, one can map the calculated (water) bead diffusion coefficients (34.5 × 10−5 cm2/s) with the experimental36 diffusion constant of water (2.43 × 10−5 cm2/s).37 This mapping provides us with a scaling factor 2 × 34.5 × 10−5/2.43 × 10−5 ≈ 28 for conversion of actual time scale to the effective time scale in the present DPD simulations. Note that since two water molecules are lumped into a single DPD bead, we have multiplied the ratio of bead to experimental diffusion coefficients by a factor of 2.38 Scaling the time unit of DPD with this time scaling factor, the relaxation time for the intermittent correlation function increases to ≈300 ns, which is comparable with the experimental27 relaxation times for micelle formation equilibrium (≈200 ns for a SDS solution of the same concentration). The relaxation time for the continuous correlation function increases to about a few nanoseconds, which can be compared with the experimental relaxation time for the reaction of monomers with micelles.35

5. CONCLUSIONS DPD is a powerful method for simulation of fluids of a variety of molecular shapes and complexities. However, there exist technical problems in dealing with electrostatic interactions between soft DPD beads. This is due to the divergence of Coulomb potential at short distances, which introduces artificial ion pairing between soft charged DPD beads. In this work, the point charges on DPD beads have been replaced with Gaussian smeared charge distributions. The potential energy of interaction between beads remains finite at short distances but recovers the Coulomb law at long distances. Screening the Gaussian charge distributions, with wider Gaussians of opposite charge, allows to split the electrostatic interaction into real- and reciprocal-space contributions. The method is easy to implement; only slight modifications in the implemented standard Ewald summations method are needed. Compared to the existing methods, it is shown that optimization of the screening parameter as well as the realand reciprocal-space cutoffs, leads to improvements in terms of convergence rate and computational speed. This method is validated against model test systems in the literature. The method has also been employed to study self-assembly in solutions of SDS in water. The CMC and the equilibrium concentration of free surfactants, in solutions with SDS concentrations varying from CMC to 20 times larger than CMC, are in close agreement with experiment. The microscopic structure of micelle and distributions of hydrophobic and hydrophilic groups of micelle and counterions at the water interface are in agreement with experiment. We have also examined the dynamics of monomer exchange between micelles and solution, in terms of the intermittent and continuous correlation functions for variation of micelle concentration with time. Two discrete relaxation processes, whose relaxation times differ by 2 orders of magnitude, are found. The relaxation times for each process are much shorter than experimentally relevant relaxation times (characterized in terms of exchange of surface molecules with solution and equilibrium between surfactants in the solution and micelles through diffusion of surfactants), but the ratio of the relaxation times is comparable with experiment. Scaling the DPD time unit with a time scaling factor (the ratio of calculated to experimental diffusion coefficient of water), however, we obtain relaxation times comparable with experiment.27,31,34

Figure 12. Decay of the continuous correlation function, Sn(t), for exchange of surfactant molecules between micelles of aggregation number n and solution. The sizes of the micelles are shown in the figure’s legend.

decays much faster. The decay rate of the function Sn(t) depends on the rate of exchange of the surface molecules with the solution. Our calculated relaxation times for micelles within the range of sizes varying from 25 to 55 (see Figure 12) vary between 0.1 to 0.25 ns. This time scale is shorter than the relaxation time for surfactant exchange, calculated in terms of the decorrelation of Cn(t), by 2 orders of magnitude and is nearly independent of size. Experimentally measured relaxation times34 for the fast relaxation process (reaction between micelles and monomers) are shorter than those for the slow relaxation process (formation/decomposition of the micelles) by two−three orders of magnitude. The distinct difference between these two relaxation times (observed in experiment and in our simulations) reveals distinct motion regimes of molecular exchange between micelles and solution, namely, dissociation/association of individual molecules from/to the micelles and establishment of equilibrium between surfactants in the solution and micelles through diffusion of surfactants in the solution. The shortness of our calculated relaxation times, compared to experiment, is expectable as the dynamics in all coarsegrained methods gets artificially accelerated.35 In this respect, the times quoted in our DPD simulations can qualitatively compare the dynamics of similar properties. Here, the ratio of our calculated relaxation times for functions Cn(t) and Sn(t), for large clusters which are close to experimentally observed micelle sizes, is about 100, which is comparable with the ratio of experimentally measured relaxation times27,34 (ranging from 100 to 1000). 4205


Article


■

(17) Chialvo, A. A.; Cummings, P. T. Simple transferable intermolecular potential for the molecular simulation of water over wide ranges of state conditions. Fluid Phase Equilib. 1998, 150−151, 73−81. (18) Jeon, J.; Lefohn, A. E.; Voth, G. A. An improved Polarflex water model. J. Chem. Phys. 2003, 118, 7504−7518. (19) Kiss, P. T.; Sega, M.; Baranyai, A. Efficient handling of Gaussian charge distributions: An application to polarizable molecular models. J. Chem. Theory Comput. 2014, 10, 5513−5519. (20) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (21) Anderson, R. L.; Bray, D. J.; Ferrante, A. S.; Noro, M. G.; Stott, I. P.; Warren, P. B. Dissipative particle dynamics: Systematic parametrization using water-octanol partition coefficients. J. Chem. Phys. 2017, 147, 094503. (22) Müller-Plathe, F. YASP: A molecular simulation package. Comput. Phys. Commun. 1993, 78, 77−94. (23) Kolafa, J.; Perram, J. W. Cutoff errors in the Ewald summation formulae for point charge systems. Mol. Simul. 1992, 9, 351−368. (24) Vila Verde, A.; Frenkel, D. Kinetics of formation of bile salt micelles from coarse-grained Langevin dynamics simulations. Soft Matter 2016, 12, 5172−5179. (25) Anachkov, S. E.; Danov, K. D.; Basheva, E. S.; Kralchevsky, P. A.; Ananthapadmanabhan, K. P. Determination of the aggregation nnumber and charge of ionic surfactant micelles from the stepwise thinning of foam films. Adv. Colloid Interface Sci. 2012, 183, 55−67. (26) Benrraou, M.; Bales, B. L.; Zana, R. Effect of the nature of the counterion on the properties of anionic surfactants. 1. CMC, ionization degree at the CMC and aggregation number of micelles of sodium, cesium, tetramethylammonium, tetraethylammonium, tetrapropylammonium, and tetrabutylammonium dodecyl sulfates. J. Phys. Chem. B 2003, 107, 13432−13440. (27) Polacek, R.; Kaatze, U. Monomer exchange kinetics, radial diffusion, and hydrocarbon chain isomerization of sodium dodecylsulfate micelles in water. J. Phys. Chem. B 2007, 111, 1625−1631. (28) Sanders, S. A.; Sammalkorpi, M.; Panagiotopoulos, A. Z. Atomistic simulations of micellization of sodium hexyl, heptyl, octyl, and nonyl sulfates. J. Phys. Chem. B 2012, 116, 2430−2437. (29) Itri, R.; Amaral, L. Q. Distance distribution function of sodium dodecyl sulfate micelles by x-ray scattering. J. Phys. Chem. 1991, 95, 423−427. (30) Jones, R. R. M.; Maldonado, R.; Szajdzinska-Pietek, E.; Kevan, L. Electron spin echo modulation of doxylstearic acid probes of the surface and internal structure of lithium dodecyl sulfate micelles: comparison with sodium dodecyl sulfate and tetramethylammonium dodecyl sulfate micelles. J. Phys. Chem. 1986, 90, 1126−1129. (31) Ravichandran, G.; Gopinath, D. Ultrasonic relaxation studies on micelle formation in aqueous solutions of some bile salts. J. Mol. Liq. 2014, 198, 122−127. (32) Rapaport, D. C. Hydrogen bonds in water. Mol. Phys. 1983, 50, 1151−1162. (33) Eslami, H.; Khanjari, N.; Müller-Plathe, F. Self-assembly mechanisms of triblock Janus particles. J. Chem. Theory Comput. 2019, 15, 1345−1354. (34) Lessner, E.; Teubner, M.; Kahlweit, M. Relaxation experiments in aqueous solutions of ionic micelles. 1. Theory and experiments on the system H20-sodium tetradecyl sulfate-NaClO. J. Phys. Chem. 1981, 85, 1529−1536. (35) Eslami, H.; Müller-Plathe, F. How thick is the interphase in an ultrathin polymer film? Coarse-grained molecular dynamics simulations of polyamide-6,6 on graphene. J. Phys. Chem. C 2013, 117, 5249−5257. (36) Partington, J. R.; Hudson, R. F.; Bagnall, K. W. Self-diffusion of aliphatic alcohols. Nature 1952, 169, 583−584. (37) Fraaije, J. G. E. M.; van Male, J.; Becherer, P.; Gracià, R. S. Calculation of diffusion coefficients through coarse-grained simulations using the automated-fragmentation-parametrization method and the recovery of Wilke-Chang statistical correlation. J. Chem. Theory Comput. 2018, 14, 479−485.

AUTHOR INFORMATION

Corresponding Author

*Email: [email protected]. ORCID

Hossein Eslami: 0000-0002-1990-0469 Florian Müller-Plathe: 0000-0002-9111-7786 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the support of this work by the Deutsche Forschungsgemeinschaft (DFG), project MU1412/ 25-2. H.E. also acknowledges the research committee of Persian Gulf University.

■

REFERENCES

(1) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett. 1992, 19, 155−160. (2) Español, P.; Warren, P. Statistical mechanics of dissipative particle dynamics. Europhys. Lett. 1995, 30, 191−196. (3) Groot, R. D. Electrostatic interactions in dissipative particle dynamics-simulation of polyelectrolytes and anionic surfactants. J. Chem. Phys. 2003, 118, 11265−11277. (4) Deserno, M.; Holm, C. How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines. J. Chem. Phys. 1998, 109, 7678−7693. (5) González-Melchor, M.; Mayoral, E.; Velázquez, M. E.; Alejandre, A. Electrostatic interactions in dissipative particle dynamics using the Ewald sums. J. Chem. Phys. 2006, 125, 224107. (6) Ewald, P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. (Berlin, Ger.) 1921, 369, 253−287. (7) Coslovich, D.; Hansen, J.-P.; Kahl, G. Ultrasoft primitive model of polyionic solutions: Structure, aggregation, and dynamics. J. Chem. Phys. 2011, 134, 244514. (8) Warren, P. B.; Vlasov, A.; Anton, L.; Masters, A. J. Screening properties of Gaussian electrolyte models, with application to dissipative particle dynamics. J. Chem. Phys. 2013, 138, 204907. (9) Warren, P. B.; Vlasov, A. Screening properties of four mesoscale smoothed charge models, with application to dissipative particle dynamics. J. Chem. Phys. 2014, 140, 084904. (10) Nam, K.; Gao, J.; York, D. M. An efficient linear-scaling Ewald method for long-range electrostatic interactions in combined QM/ MM calculations. J. Chem. Theory Comput. 2005, 1, 2−13. (11) Prinsen, P.; Warren, P. B.; Michels, M. A. J. Mesoscale simulations of surfactant dissolution and mesophase formation. Phys. Rev. Lett. 2002, 89, 148302. (12) Lee, M.-T.; Vishnyakov, A.; Neimark, A. V. Calculations of critical micelle concentration by dissipative particle dynamics simulations: The role of chain rigidity. J. Phys. Chem. B 2013, 117, 10304−10310. (13) Johnston, M. A.; Swope, W. C.; Jordan, K. E.; Warren, P. B.; Noro, M. G.; Bray, D. J.; Anderson, R. L. Toward a standard protocol for micelle simulation. J. Phys. Chem. B 2016, 120, 6337−6351. (14) Eslami, H.; Bahri, K.; Müller-Plathe, F. Solid-liquid and solidsolid phase diagrams of self-assembled triblock Janus nanoparticles from solution. J. Phys. Chem. C 2018, 122, 9235−9244. (15) Mao, R.; Lee, M.-T.; Vishnyakov, A.; Neimark, A. V. Modeling aggregation of ionic surfactants using a smeared charge approximation in dissipative particle dynamics simulations. J. Phys. Chem. B 2015, 119, 11673−11683. (16) Anderson, R. L.; Bray, D. J.; Del Regno, A.; Seaton, M. A.; Ferrante, A. S.; Warren, P. B. Micelle formation in alkyl sulfate surfactants using dissipative particle dynamics. J. Chem. Theory Comput. 2018, 14, 2633−2643. 4206


Article

Journal of Chemical Theory and Computation (38) Groot, R. D.; Rabone, K. L. Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophys. J. 2001, 81, 725−736.

4207


Gaussian charge distributions for incorporation of electrostatic

Recommend Documents