Available online at www.sciencedirect.com
Physica A 319 (2003) 421 – 431
www.elsevier.com/locate/physa
Monte Carlo simulations for a model of amphiphiles aggregation M. Girardi∗ , W. Figueiredo Departamento de F sica, Universidade Federal de Santa Catarina, 88040 900, Florian opolis, Santa Catarina, Brazil Received 27 May 2002
Abstract In this work we employed Monte Carlo simulations to study a spin-like lattice model that resembles some properties observed in micellar systems. This simple model for molecular aggregation was studied in dimensions d = 1; 2 and 3. The interaction energy favors the formation of a well-de0ned optimal aggregate of size n, where n is the number of possible states of the amphiphiles. In this model, both the amphiphiles and the water molecules (we treated aqueous solutions) occupy only a single site in the lattice. The simulations and numerical calculations in one dimension agree with the exact ones obtained in previous works. The minimum and the maximum in the aggregate-size distribution curve and the critical micellar concentration, which are the 0ngerprints of micellar aggregation, are obtained. The two- and three-dimensional cases also exhibit these features. We calculated the exponent that measures the dependence of the di4erence in height between the maximum and the minimum in the aggregate-size distribution curve on temperature. We found ∼ = 1, for all values of n we considered, independently of d. c 2002 Elsevier Science B.V. All rights reserved. PACS: 64:75: + g; 82.60.Lf; 64.60.ht Keywords: Micelles; Critical micellar concentration; Monte Carlo simulations
1. Introduction The study of micellar aggregation is of great interest not only for the industrial purposes but also to the biological sciences and medicine [1–7]. In the last years, many e4orts were made to obtain information about the size, shape and the formation ∗
Corresponding author. E-mail addresses:
[email protected] (M. Girardi),
[email protected] (W. Figueiredo).
c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/03/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 1 3 9 0 - 0
422
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
mechanism of these structures. The use of extensive computer simulations played an important role in this task [8–13]. In previous studies, we considered a simple chain model to describe amphiphilic aggregation in two and three dimensions through Monte Carlo simulations [14,15]. The amphiphilic molecules were mapped as chains of connected sites. One site, with hydrophilic properties, represents the surfactant head, and the other sites model the hydrophobic part of the molecule chain. For this model we found the critical micellar concentration (CMC) and de0ned the exponent associated with the vanishing of the height di4erence between the maximum and the minimum in the aggregate-size distribution curve (ADC). This exponent measures of the stability of the most probable aggregate against temperature variations. From the simulations we found ∼ = 1:0 (d = 2) and ∼ = 2:0 (d = 3). It was also noted that, even for a total concentration of amphiphiles above the CMC, the ADC can be a monotonically decreasing function of the aggregate size, for temperatures above a micellization temperature. We also performed simulations and obtained analytical results for a spin-1 lattice model [16,17] where the amphiphiles are mapped onto its ±1 spin components, and the water molecules are represented by its 0 spin component. This model does not exhibit the minimum and the maximum in the aggregation curve [18] and cannot be used to mimic the micellar aggregation. However, it was shown that, choosing a di4erent set of interaction energies [19], the system behaves like a micellar solution where the both extrema in the ADC are present. In this paper, we extend the simple one-dimensional model of molecular aggregation [19] to two and three dimensions, employing Monte Carlo simulations. The amphiphilic molecules occupy a single site in the lattice and possess a discrete set of internal degrees of freedom. The molecular interactions favor the formation of clusters with a well-de0ned optimal size. This simple model incorporates the main properties observed in the amphiphilic solutions in equilibrium. A local minimum and maximum in the aggregate-size distribution curve and a critical micellar concentration for suitable values of temperature and total concentration of amphiphiles. An equivalent model was exactly solved by Duque and Tarazona [20] in a continuous medium, where the amphiphiles were represented by rods. In their model, the energy couplings depends on the distance between two neighboring rods and of their respective states (related to their internal degrees of freedom of the rod). This work is organized as follows: in Section 2 we describe the model, and a suitable energy function is proposed to describe it in one dimension. In Section 3 we exhibit some exact results in one dimension obtained through the transfer matrix formalism. For the two- and three-dimensional versions of the model we present only simulation results. In Section 4 we discuss our results and we draw our conclusions. 2. The model Let us initially consider the one-dimensional version of the model. In this lattice model, both the amphiphilic and solvent molecules occupy a single site in the lattice, instead of the usually studied chains model, where the surfactant molecules occupy many sites. The amphiphiles possess internal degrees of freedom with n di4erent states.
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
423
These states can be associated to the various conformations and sizes of the amphiphiles in the real systems. We expect that increasing n the resulting aggregates become large as seen for the micelles when the number of monomers of the amphiphiles grows. For the water molecules only a single state is assigned. For the one-dimensional case, we can write the following energy function for the state {qi }: N {−J [ (1+qi ;qi+1 ) + (qi ;n) (qi+1 ;0) ] − (qi ;0) } ; (1) H= i=1
where qi is one of the n possible states (qi = 1; 2; : : : ; n) of the amphiphile in the ith site or 0 if there is a water molecule in the ith site. functions are the delta’s Kronecker, N is the lattice size, J is the coupling constant (we assume J ¿ 0) and is the chemical potential of the water. By changing we can control the degree of dilution of the system. From Eq. (1), it is clear that sequences like 012 : : : n0 are energetically favored, and we may expect that the ground state of the model is of the type: : : : 012 : : : n012 : : : n0 : : : . It means that, for very low temperatures and = 0, the system presents only amphiphilic aggregates of size n, separated by one water molecule. Otherwise, for any other choice of the ordering states of the molecules at zero temperature, we would have aggregates of in0nite size as a large droplet of oil inside water. This trend to form aggregates with some well-de0ned size is the main property of this model. Micellar aggregates are de0ned by a set of neighboring amphiphiles, in any state, surrounded by water molecules. The total concentration of amphiphiles can be obtained by the equation 1 Xt = (1 − (qi ;0) ) ; (2) N i where the sum is over all the sites and N is the total number of lattice sites. The concentration of aggregates of size m is given by Nm Xm = ; (3) N where Nm is the number of aggregates of size m. The energy function given by Eq. (1) is not a suitable form to represent the energies in two and three dimensions. Thus, we have to employ some tricks to de0ne the aggregates in higher dimensions. Fig. 1 shows a typical aggregate of size m = 5 on a square lattice with the amphiphiles labelled by their states. Let us take this con0guration as being that of minimal energy. With this choice we are keeping the same features as presented in the one-dimensional case: the stacking repeating blocks, at very low temperatures, have n particles, where n is the number of the possible states of the amphiphiles. In this way, we can represent the interaction energies by the matrices of Table 1, one for each axis. One can read these energies from the matrix according to the following rules: in the X (Y ) axis matrix, the 0rst column indicates particles that are in the x(y) position while the 0rst line indicates particles in x + 1(y + 1) position. As stressed above, the aggregates of minimal energy must contain n amphiphiles, each one in a di4erent state. Otherwise, aggregates of in0nite size would appear at zero temperature. Di4erent forms and sizes of aggregates are obtained by changing the number n of states and the values of the couplings in the matrices of Table 1. This
424
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
0
0
0
0
0
0
0
4
0
0
0
1
2
3
0
0
0
5
0
0
0
0
0
0
0
Fig. 1. Low temperature most probable aggregate for n = 5 in two dimensions.
Table 1 Matrix representation for the couplings in X and Y directions for n = 5 in two dimensions X -axis
0 1 2 3 4 5
Y -axis
1
2
3
4
5
0
−J 0 0 0 0 0
0 −J 0 0 0 0
0 0 −J 0 0 0
−J 0 0 0 0 0
−J 0 0 0 0 0
0 0 0 −J −J −J
0 1 2 3 4 5
1
2
3
4
5
0
−J 0 0 0 0 0
0 0 0 0 −J 0
−J 0 0 0 0 0
−J 0 0 0 0 0
0 0 −J 0 0 0
0 −J 0 −J 0 −J
choice seems arbitrary, because for a given n and coupling constant the aggregate, at the ground state, is unique. However, for non-zero temperatures the most probable aggregates exhibit di4erent sizes and shapes. As we can see in the next section, these aggregates present a number of molecules less than n for 0nite temperatures. For the three-dimensional case we simply add another matrix to be associated to the z-axis. 3. Exact results and simulations Let us initially consider the one-dimensional version of the model presented in the last section. Using the transfer matrix technique [21], we can calculate the partition function of the model and determine, for each value of temperature and chemical potential, the total concentration of amphiphiles and the aggregation distribution curve. The transfer matrix in the case of n states is a non-symmetric (n + 1) × (n + 1) matrix.
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
425
Although some eigenvalues are complex, the pair of complex conjugate eigenvalues N and (? )N render a positive contribution to the partition function when they are added. A similar calculation can be seen in Ref. [16] and we will not present the details here. However, we observed that only for = 0, a closed expression for the eigenvalues is available. With this value of the chemical potential, the average concentration of amphiphiles is exactly n=(n+1). We are interested in the regime of low dilution (where the micellar aggregates occur), which is not the case for the value = 0. For instance, it is easy to see that, for n=1, Eq. (1) is the proper Ising Hamiltonian, and the fraction of amphiphiles (indexed by up spins) is one half. In order to access the regime of low concentrations we need to consider values of di4erent from zero. Unfortunately, for these values of , the expressions are not easily manageable and we need to resort to the numerical calculations. We have also employed Monte Carlo simulations to investigate this one-dimensional model. The simulations were carried out in a lattice with 104 sites and we used periodic boundary conditions. As a 0rst step in the simulation, for a given total concentration of amphiphiles Xt , we distributed them randomly in the lattice, and their states were also chosen randomly as an integer number between 1 and n. The remaining sites of the lattice were set to 0, indicating the presence of water molecules. The initial state was chosen to be completely random. Two di4erent dynamics were used to describe the time evolution of the system throughout the phase space: one is the Kawasaki exchange dynamics [22] where we choose, at random, two di4erent sites and exchange the particles sitting on them (water or amphiphile) according to the Metropolis [23] prescription. This dynamical process mimics the di4usion of the particles throughout the lattice. The other dynamical process is the single spin-Kip Glauber kinetics [24] where an amphiphile is chosen at random in the lattice and its state can change to another one, between 1 and n, accordingly to the Metropolis rule. We must emphasize that the total concentration is never modi0ed during the simulation. We de0ne one Monte Carlo step (MCs) as been N = Ld trials of executing both dynamics, where d is the spatial dimensionality and L is linear dimension of the lattice. To calculate the thermodynamic averages, the system must 0rst reach the equilibrium, which occurs typically after 5 × 104 MCs in one dimension. The averages are then calculated performing more 105 MCs. In Fig. 2 we exhibit the concentration of the free amphiphiles as a function of the total concentration in one, two and three dimensions. In one dimension, we show the exact numerical results and the simulations performed by employing the two dynamical processes with the same weight. These are represented by the continuous line (exact result) and the full circles (simulations). As to be expected, this curve follows an ideal gas law at very low concentrations. However, at higher concentrations, the amphiphiles prefer to aggregate and the curve departs from the ideal behavior. Also, this signals that a CMC was already attained. Fig. 2 shows the results of simulations for two (connected triangles) and three (connected crosses) dimensions. As we can see in this 0gure, as the total concentration increases, the aggregation is easier in d = 3, that is, the CMC is attained at lower concentrations. Fig. 3 displays the numerical results obtained for the di4erence in height between the maximum and the minimum in the aggregate-size distribution curve as a function of temperature in one dimension. These values, plotted near the temperature TM = 0:3355
426
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
X1
0.2
0.1
0 0
0.1
0.2
Xt Fig. 2. Concentration of isolated amphiphiles X1 as a function of the total concentration Xt . Exact numerical calculations (continuous line) and simulations (circles) for N = 104 ; n = 4 and T = 0:3 in the one-dimensional case. Simulation in two dimensions (connected triangles) for n = 5 and T = 0:256, and in three dimensions for n = 8 (cubic shape) and T = 0:400 (connected crosses). The dashed line gives the ideal gas behavior.
(the temperature is normalized in units of J=kB ), are well 0tted to a straight line, that is, ∼ (TM − T ) , where = 1:0 and TM is the micellization temperature. The value of we found from the simulations for this unidimensional model is =1:0±0:1. The ADC shown in the inset of Fig. 3, exhibits the minimum and the maximum characteristic of the micellar formation for temperatures below TM . Above this temperature the ADC is a monotonically decreasing function of the aggregate size. Figs. 2 and 3 were also built for other values of n, besides n = 4. For all these values of n we have found the same value for the exponent . The agreement between numerical calculations and Monte Carlo simulations in one dimension, prompted us to perform Monte Carlo simulations on this model in higher dimensions. It is very easy to implement the two- and three-dimensional simulations for this model. Unfortunately, we were not able to write a suitable energy function for dimensions di4erent from 1. For the set of interaction energies shown in Table 1, the ground state of the model is well-de0ned. However, a slightly change in the dynamical process that drives the system to equilibrium is necessary. We do not perform the Glauber dynamics when the aggregates are non-linear. In this way, we keep 0xed the number of amphiphiles with the same state during the simulations. At the beginning of the simulation, we select at random the positions of all amphiphiles for a given total concentration. Besides, the initial distribution of states is chosen to be uniform. The use of the Glauber kinetics to the formation of two- and three-dimensional aggregates is not so eLcient as in the case of one-dimensional aggregates. This happens because the
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431 0.005
427
0.6
0.4 mX m
0.004
0.2 ∆
0.003
∆
0 0
2
0.002
4 m
6
8
0.001
0 0.33
0.332
0.334
0.336
T Fig. 3. Plot of the parameter as a function of temperature T in one dimension for n = 4 and Xt = 0:3. Exact results (dashed line) and simulations (circles) clearly show a linear behavior. The inset shows the aggregate-size distribution curve for two di4erent temperatures, T =0:300 (crosses) and TM =0:3355 (circles). It also shows the di4erence in height for the curve with T = 0:300. At TM the di4erence is null.
energetics we impose to the formation of general aggregate shapes in higher dimensions turn the relaxation times too long. In one dimension this is not the case. Then, in order to speed up the simulations in two and three dimensions we have considered only the Kawasaki dynamics, where the number of amphiphiles in each state is conserved during all the steps of the simulation. The Monte Carlo simulations were carried out for lattices with linear sizes L = 100 and 20 in two and three dimensions, respectively. We present the results in two and three dimensions and for three di4erent number of states: n=5; 8 and 12. In two and three dimensions we can built, for a given n, aggregates of di4erent shapes. In this work we take for n = 5 only Kat aggregates, for n=8, Kat and cubic shape aggregates, and for n=12, Kat and ellipsoidal aggregates. For d = 2 and 3, the thermalization process is slower than in one dimension, and we need 104 MCs to reach the equilibrium state. Fig. 4 presents the ADC in d = 2 and 3 for two di4erent temperatures. Clearly, we observe that for higher temperatures the curves are monotonically decreasing, while for lower values of T , the ADC presents a local minimum and maximum characteristic of the micellar systems. As pointed out before, the largest aggregate sizes are obtained only at very low temperatures and it coincides with the number of states n of the amphiphiles. For instance, in Fig. 4(a), where d = 2; n = 5 and T = 0:256 (crosses), the optimal aggregate size is 3. Fig. 5 displays the plots of the di4erence in height between the maximum and the minimum in the ADC as a function of temperature for three di4erent values of n, for the two-dimensional case. In these plots the curves are 0tted to ∼ (TM − T )
428
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
0.8
(a)
(b)
mX m
0.6
0.4
0.2
0 2
0
4
8 0
6
2
4
6
8
10
m Fig. 4. Aggregate-size distribution curves for two di4erent temperatures. (a) d = 2; n = 5 and Xt = 0:1. Circles (T = 0:400), crosses (T = 0:256). (b) d = 3; n = 8 (cubic shape) and Xt = 0:1. Circles (T = 0:500), crosses (T = 0:380). The lines serve to guide the eyes.
0.006
(a)
(c)
(b)
∆
0.004
0.002
0 0.300
0.304
0.300 0.303
0.306 0.309
T Fig. 5. Plot of the parameter as a function of temperature T in two dimensions for n = 5 (a), n = 8 (b) and n = 12 (c). The concentration is Xt = 0:1. Simulations are given by the circles and the straight lines are M (b) = 1:0M and (c) = 1:1. M the best 0t to the data points. The values of the exponent are: (a) = 1:1,
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
(a)
0.01
(b)
429
0.0025
0.008
0.002
0.006
0.0015
0.004
0.001
0.002
0.0005 0
0
∆
0.3 0.302 0.304 0.306
0.435
(c)
0.006
0.438
0.441
(d) 0.003
0.004
0.002
0.002
0.001 0
0 0.3
0.303 0.306
0.308 0.31 0.312 0.314
T Fig. 6. Plot of the parameter as a function of temperature T in three dimensions for (a) n = 5, (b) n = 8 (aggregate of cubic shape), (c) n = 8 (Kat shape) and (d) n = 12 (ellipsoidal shape). The concentration is Xt = 0:1. Simulations are given by the circles and the straight lines are the best 0t to the data points. The M exponent is the same for all plots, = 1:1.
where TM is the micellization temperature. Above this temperature TM , the system does not show micellar aggregates. The values of are near 1 for all the values of n we investigated. Fig. 6 shows also as a function of temperature for three di4erent values of n in three dimensions. Fig. 6(a) is for n = 5, Fig. 6(b) n = 8 (cubic aggregate), Fig. 6(c) n = 8 (Kat aggregate) and Fig. 6(d) n = 12 (ellipsoidal aggregate). We see that, independently of n and of the shape of the aggregate, is near 1. From the exact results in one dimension, and the simulations performed in two and three dimensions, we observed that the exponent does not depend neither on dimensionality nor on the number of con0gurations of the amphiphiles. Finally, in Fig. 7, we exhibit the behavior of the parameter as a function of the reduced temperature (|T − TM |=TM ) for n = 8 in one, two and three dimensions for a 0xed concentration, and for temperatures well below the micellization temperature. Far from TM , as expected, the aggregation in three dimensions is enhanced relatively to the aggregation in one and two dimensions. 4. Discussions and conclusions In real amphiphilic systems, the carbon chain size and the hydrophilic properties of the surfactant heads, play the most important role on the size and shape of the micelles. For example, strongly ionic heads induce the formation of spherical micelles due to its large headgroup area [25]. In this work, we assembled all the relevant microscopic information of the amphiphiles in the parameter n, which mimics all the
430
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431 0
10
−1
∆
10
−2
10
−3
10
−4
10
0
0.1
0.2
0.3
0.4
|T-TM | TM Fig. 7. Logarithmic plot of the parameter as a function of the reduced temperature |T − TM |=TM in one (circles), two (crosses) and three (triangles) dimensions. The concentration is the same for all dimensions (Xt = 0:3) and n = 8. For the three-dimensional case the shape of the aggregate is cubic. The lines serve to guide the eyes.
conformational states of the surfactants, and in the coupling constants of the interaction matrices. This simpli0cation saves the simulation CPU time, and, at same time, renders a much simpler algorithm than the usual one employing the reptation movements. We have employed Monte Carlo simulations to 0nd the aggregate-size distribution curve for the model in one-, two- and three dimensions. In the particular case of the one-dimensional version of this model, transfer matrix techniques were also used and some exact results were determined. The results we obtained for the concentration of free amphiphiles as a function of the total concentration, and for the aggregate-size distribution curve capture the main properties exhibited by the real micellar systems. We have calculated the exponent , which is related to the vanishing of the di4erence in height between the maximum and the minimum in the ADC, as a function of temperature. This exponent is not related to any critical phenomenon, but its importance here is to account for the stability of the most probable aggregates. Larger the exponent more stable aggregates are found. We have determined the value ∼ = 1 for all the investigated values of n and for dimensions d = 1; 2 and 3. This result is not the same as observed in the chain models of aggregation of amphiphiles. For instance, our Monte Carlo simulations [14,15] for chain models in two and three dimensions give = 1 and 2, respectively. As the packing of the chains is more e4ective in three than in two dimensions, then the exponent must be larger in three than in two dimensions.
M. Girardi, W. Figueiredo / Physica A 319 (2003) 421 – 431
431
On the other hand, in the present model, the exponent is insensitive to both spatial dimension and number of possible conformational states of a surfactant. We think this behavior is due to the fact that, to change the internal energy of an aggregate, in any dimension, the most probable path is through the addition or withdrawn of a single particle. Indeed, in real chain amphiphiles we can easily change the internal energy of an aggregate without changing the number of its particles. Acknowledgements This work was partially supported by the Brazilian agency CNPq. References [1] C. Tanford, Hydrophobic E4ect: Formation of Micelles and Biological Membranes, Wiley, New York, 1980. [2] W. Gelbart, A. Ben-Shaul, D. Roux (Eds.), Micelles, Membranes, Microemulsions, and Monolayers, Springer, New York, 1994. [3] Y. Chevalier, T. Zemb, Rep. Prog. Phys. 53 (1990) 279. [4] A.Z. Panagiotopoulos, M.A. Floriano, S.K. Kumar, Langmuir 18 (2002) 2940. [5] U. Reimer, M. Wahab, P. Schiller, H.J. MRogel, Langmuir 17 (2001) 8444. [6] M. LSTsal, C.K. Hall, K.E. Gubbins, A.Z. Panagiotopoulos, J. Chem. Phys. 116 (2002) 1171. [7] G. Cristobal, J. Rouch, J. CurSely, P. Panizza, Physica A 268 (1999) 50. [8] P.G. Bolhuis, D. Frenkel, Physica A 244 (1997) 45. [9] R.G. Larson, J. Phys. II France 6 (1996) 1441. [10] D. Brindle, C.M. Care, J. Chem. Soc. Faraday Trans. 88 (1992) 2163. [11] J.C. Desplat, C.M. Care, Mol. Phys. 87 (1996) 441. [12] P.H. Nelson, G.C. Rutledge, T.A. Hatton, J. Chem. Phys. 107 (1997) 10 777. [13] A.T. Bernardes, V.B. Henriques, P.M. Bisch, J. Chem. Phys. 101 (1994) 645. [14] J.N.B. de Moraes, W. Figueiredo, J. Chem. Phys. 110 (1999) 2264. [15] M. Girardi, W. Figueiredo, J. Chem. Phys. 112 (2000) 4833. [16] M. Girardi, W. Figueiredo, Phys. Rev. E 62 (2000) 8344. [17] M. Girardi, W. Figueiredo, Phys. Stat. Sol. A 187 (2001) 195. [18] H. Wennerstrom, B. Lindman, Phys. Rep. 52 (1979) 1. [19] D. Duque, Phys. Rev. E 64 (2001) 063 601. [20] D. Duque, P. Tarazona, J. Chem. Phys. 107 (1997) 10 207. [21] C.J. Thompson, Mathematical Statistical Mechanics, Princeton University Press, Princeton, NJ, 1972. [22] K. Kawasaki, in: C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, Vol. 4, Academic, London, 1974. [23] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [24] R.J. Glauber, J. Math. Phys. 4 (1963) 294. [25] J. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1992.
