Available online at www.sciencedirect.com

Physica A 324 (2003) 621 – 633

www.elsevier.com/locate/physa

Lattice-gas with two- and three-body couplings as a model for amphiphilic aggregation M. Girardi∗ , W. Figueiredo Departamento de F sica, Universidade Federal de Santa Catarina, 88040-900, Florian opolis, Santa Catarina, Brazil Received 17 October 2002

Abstract In this paper we studied a lattice-gas model with two- and three-body interactions in linear, square and cubic lattices. We observed that, at low concentrations, the model presents some aggregation characteristics similar to micellar solutions. These peculiarities include a plateau in the plane of concentration of isolated molecules as a function of the total concentration (CMC), and a local minimum and maximum in the aggregate-size distribution curve (ADC). The transfer matrix technique, Monte Carlo simulations, and the independent cluster approximation were employed to 4nd the ADC, CMC and the micellization temperature of the model. The one-dimensional case was solved exactly and a transition between a micellized and a non-micellized state is displayed. For the two- and three-dimensional versions of the model, extensive Monte Carlo simulations were performed in order to 4nd the exponent
associated with the di5erence between the local minimum and maximum heights in the aggregate-size distribution curve. We have found that
is 1, independent of the spatial dimension. c 2003 Elsevier Science B.V. All rights reserved.  PACS: 64.75.+g; 82.60.Lf; 64.60.ht Keywords: Micelles; Independent cluster approximation; Monte Carlo simulations

1. Introduction Amphiphilic molecules consist, in general, of a hydrophobic carbonic chain connected to a hydrophilic head group. These molecules have special properties that allow to an increase of the miscibility of oil in water [1–3], and favour their aggregation ∗

Corresponding author. E-mail addresses: [email protected] (M. Girardi), [email protected] (W. Figueiredo).

c 2003 Elsevier Science B.V. All rights reserved. 0378-4371/03/$ - see front matter
doi:10.1016/S0378-4371(03)00079-7

622

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

in a broad variety of structures like micelles, vesicles and membranes. Due to these facts, a large e5ort was made in the last years to understand the physics behind the rich behaviour of the surfactants in solution. Various lattice and o5-lattice models for the aggregation of amphiphiles in aqueous solutions were investigated to 4nd the equilibrium properties of the solutions at low and high amphiphilic concentrations. In particular, at low concentrations, the appearance of a special type of aggregate, called micelle, was largely investigated [4–18]. Micellar systems are characterized by the presence of clusters of amphiphiles that are in chemical equilibrium with isolated amphiphiles in the solution. In the thermodynamical equilibrium, there is a continuous exchange of amphiphilic molecules between clusters of di5erent sizes. A look at the distribution curve of the aggregate-sizes in a micellized solution, reveals the presence of a maximum in the concentration of the isolated amphiphiles, and another one at a preferred size [18]. This behaviour of the aggregate-size distribution curve and the plateau observed in the curve of the free amphiphile concentration versus total concentration are the 4ngerprints of the micellar phase. Various theoretical approaches based on mean 4eld calculations [19] and powerful computer simulations [4–14] were employed to shed some light into the mechanisms of the micellar formation, mainly focusing the size, shape and stability of the micelles. However, the complexity of the system prevent us to extract more detailed results. This complexity arises because the amphiphilic molecule presents a polymer like structure, where a large amount of monomers, some of them of the hydrophobic type, and some other of hydrophilic type, are connected forming a Fexible chain. It means that, a theoretical treatment of such amphiphilic system needs to take into account the conformational degeneracy of the molecules, which is not a simple task. Even the use of simulations becomes a hard job. Molecular dynamics requires a detailed modelling of the interaction potential among the monomers and we must solve a large number of equations of movement to obtain the properties of the system. Monte Carlo simulations are simpler to implement, but we need long time runs and lots of samples are necessary to reduce the statistical errors. Therefore, in the present work, we propose a simple lattice-gas model that exhibits some important properties of the micellar solutions. Lattice-gas models are easier to treat theoretically, since the molecules have no structure. The hamiltonian of the present model includes two- and three-body interactions and the system is examined in the thermodynamical equilibrium. This lattice-gas model was studied in the last years, but its aggregation characteristics at low concentrations were not explored at all. Chin and Landau [20], for example, studied the triangular Ising lattice-gas with two- and three-body interactions, employing Monte Carlo simulations. They obtained the phase diagram in the 4eld-temperature and the coverage-temperature planes for a wide range of the interaction parameters. Some exact results were also obtained by Baxter and Wu [21] for the Ising model (that can be mapped onto the lattice-gas model by a simple change of variables) in the triangular lattice with three-body interactions, and by Wu and Wu [22] in the KagomIe lattice with two- and three-site couplings. However, no information concerning the aggregation behaviour at low concentrations was considered. In this way, the aim of the present paper is to extract some relevant data from this lattice-gas model that can help us to understand the mechanisms of the

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

623

micellar aggregation. In our model, the amphiphiles have no internal degrees of freedom, occupying a single site in the lattice. The other unoccupied sites are associated with water molecules. These over simpli4cations allow us to obtain some exact results in one dimension by the transfer matrix technique, and to perform accurate Monte Carlo simulations for the two- and three-dimensional versions of the model. We also consider an independent cluster approximation to compare with the results of the simulations. This approximation was originally employed to investigate the properties of the imperfect gases [23]. We calculated the aggregate-size distribution curve, the critical micellar concentration and the exponent
, which measures the stability of the aggregates under temperature variations. This paper is organized as follows: in Section 2 we describe the model and de4ne the quantities of interest. In Section 3 we exhibit the exact calculations in one dimension, where we apply the transfer matrix technique, and we also introduce the independent cluster approximation. In Section 4 we describe the Monte Carlo simulations, and 4nally, in Section 5, we present our results and conclusions. 2. The model In this paper we studied a lattice-gas model with two- and three-body interactions. In this model, each particle of the gas represents one amphiphilic molecule, while the holes are associated with the water molecules of the solution. It means that neither the amphiphilic molecules nor the water molecules have any internal degrees of freedom or spatial structure. We treated here only dilute solutions for which the total concentration of the gas (Xt ) never exceeds 0:1, in order to observe the formation of micelles. The hamiltonian of the system is given by the equation

H = −J ni nj − D ni nj nk ; (1) {i; j}

{i; j; k}

where the 4rst sum is taken over all pairs of nearest-neighbors, and the second one is over all triplets of neighboring molecules disposed in a line [24]. J and D are the coupling constants and ni = 0; 1 denotes the occupation variable: ni = 0 represents a water molecule and ni = 1 an amphiphilic molecule. We call attention to the fact that, according to Eq. (1), the water–water and water– amphiphile interactions, contribute equally to the energy of the system, and are chosen to be zero, because we assigned the variable ni = 0 to the water molecules. If we had chosen another value for the interactions of the water molecules, the only e5ect would be a shift in the value of the coupling constant J , and the addition of a chemical potential for the amphiphiles. The two- and three-body interactions in Eq. (1) aim to represent the surface and volume energies of the micellar aggregates, which, in a real system, compete to limit the size of aggregates. As discussed in Refs. [2,3], the formation of micellar aggregates in water is related to a delicate balance of the following mechanisms: (a) the hydrophobic effect, which segregates the carbonic chains of the amphiphiles, (b) the repulsion among the head-groups of the amphiphiles, which limits the size and shape of the formed

624

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

micelles, and (c) the aKnity of the head-groups with the water. The simple model given by the hamiltonian of Eq. (1), is able to reproduce the main characteristics of the micellar aggregates by choosing the parameters J and D to have suitable values. To demonstrate the complex behaviour that arises from the hamiltonian, Eq. (1), we consider the ground-state of the system, when the temperature T is zero. For J ¿ 0 and D ¿ 0, the aggregation of the gas particles is favoured, and a large aggregate containing all amphiphiles appears. The same happens when J ¿ 0 and J |D| or D ¿ 0 with D|J |. On the other hand, for J ¡ 0 and D ¡ 0 (at low concentrations), the particles like to be completely dispersed in the solution. However, out of this range of values of J and D, the behavior of the system is not trivially determined, even at T = 0. For instance, in one dimension, when |D|2J with J ¿ 0 and D ¡ 0, and for a concentration Xt ¡ 23 , the molecules aggregate only in pairs, and the ground-state is degenerate, containing all possible states in which the molecules form pairs. Exactly at Xt = 23 , the ground-state is unique, and we have pairs of molecules separated by only a single hole. We can write the total concentration of amphiphiles as follows: V ni N Xt = i=1 = ; (2) V V where N is the total number of amphiphiles and V is the total number of sites de4ned by V = Ld , L is the linear lattice size, and d is the dimensionality of the lattice. The concentration of aggregates of size m is written as Xm =

Nm ; V

(3)

where Nm is the number of aggregates of size m. We de4ne an aggregate of amphiphiles of size m as a set of m nearest-neighbor sites with ni =1. Another quantity of interest is the exponent
that measures the thermal stability of the micellar aggregates. The curve mXm versus m, called aggregate-size distribution curve (ADC), displays a maximum and a minimum when the amphiphilic system is micellized [18]. As the temperature of the system increases, the di5erence in height  between the maximum and the minimum in the ADC vanishes, and this curve becomes monotonically decreasing with the aggregate size. We associate a power law to this vanishing in the form  ˙ (T − TM )
, where we de4ned TM as been the micellization temperature. Above this temperature, we have no more micelles. We would like to stress that we have no proper critical phenomenon here, and the exponent
gives only a measure of the stability, under thermal Fuctuations, of the micellar aggregates.

3. Calculations In the 4rst step, we derive some exact results for the simpler case where the amphiphiles are con4ned in a linear lattice with periodic boundary conditions. This one-dimensional version of the model can be solved exactly by the transfer matrix technique [25]. From the hamiltonian of Eq. (1), we can write the grand-partition

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

625

function for the model =

V





exp(N )

N =0 V


=

exp(−H )

{ni }

R(ni−1 ; ni | nˆi nˆi+1 ) ;

(4)

{ni ;nˆi } i=1

where  is the chemical potential of the amphiphiles,  =(kB T )−1 , kB is the Boltzmann constant and R(ni−1 ; ni | nˆi nˆi+1 ) is an element of the 4 × 4 transfer matrix R, which is given by   J J ni−1 ni + nˆi nˆi+1 R(ni−1 ; ni | nˆi nˆi+1 ) = ni ;nˆi exp 2 2 ×exp(Dni−1 ni nˆi nˆi+1 + ni nˆi ) :

(5)

In the last equation we used the identity ni = n2i and applied periodic boundary conditions. The new de4ned variable nˆi = 0; 1 becomes equal to ni due to the delta’s Kronecker ni ;nˆi . To obtain some information about the degree of aggregation of the system, we must write the probability to have a cluster of size m. This cluster of size m is de4ned by a sequence of m amphiphiles between two water molecules (0; 1; 1; : : : ; 1; 1; 0). This probability is written as 1 P(m) = (1 − ni ;1 )ni+1 ;1 · · · ni+m+1 ;1 (1 − ni+m+2 ;1 ) ; (6)  where we have ni ;1 = ni and the meaning of the brackets is  =

V


R(ni−1 ; ni | nˆi nˆi+1 ) :

(7)

{ni ;nˆi } i=1

The occupation variable must be written as a 4×4 matrix. This is done by the tensorial product of the 2 × 2 identity matrix and the usual occupation matrix as   1 0 0 0       0 0 0 0 1 0 1 0   ni ≡ n = ⊗ = (8)  : 0 0 1 0 0 1 0 0   0 0 0 0 We can rede4ne Eq. (6) using the matrix notation, and it becomes 1 P(m) = Tr[(U −1 RU )L−m−1 U −1 (1 − n)Rn : : : nR(1 − n)U ] ; (9)  where U is the matrix that diagonalizes the transfer matrix R and 1 is the identity matrix. In the thermodynamic limit V → ∞, Eq. (9) gives the concentration of aggregates of size m. The total concentration of amphiphiles as a function of the chemical potential  can be easily obtained di5erentiating Eq. (4) relative to . We will not show here the explicit expressions for P(m) and Xt since they are very lengthy.

626

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

2J+D

2J

2J 2J

2J+D

2J

Fig. 1. The six possible arrangements of three molecules that form an aggregate of size m = 3 with their respective energies as a function of the parameters J and D.

The determination of the matrix U and all the other mathematical manipulations were done with the aid of the symbolic calculation software Maple. The exact treatment of the hamiltonian given in Eq. (1) in two and three dimensions is not possible as we did in one dimension, and we must resort to another approach. In this work we apply the formalism used in the context of the imperfect gases, which is based on the independent cluster approximation. In this approximation it is assumed that clusters of various sizes exist and that there is a chemical equilibrium between these clusters. No cluster–cluster interaction is considered. De4ning the partition function and the chemical potential for a cluster of size m as qm and m , respectively, and considering the chemical equilibrium condition m =m ( is the chemical potential of a single molecule), we 4nd the equilibrium aggregate-size distribution function Xm = qm exp(m) :

(10)

The partition function qm is obtained by counting the number of arrangements of m molecules to form a cluster of size m, where each molecule has at least one bound, multiplied by the Boltzmann factor exp(E) (E is the energy of the cluster). For instance, Fig. 1 shows the six possible clusters for m = 3 in a square lattice, and their energies. Our counting of the clusters is similar to that employed in the lattice animals studies, and it can be easily implemented applying the algorithm proposed by Martin [26] or an improved version given by Mertens [27,28]. In one dimension, we have only one way to arrange m molecules to form a cluster; thus, the partition function qm

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

627

is straightforward and we can write the expression for Xm as X1 = exp() ; Xm = m exp[(m − 1)J + (m − 2)D + m];

(11) m¿2

(12)

and the total concentration is Xt =

∞


mXm = {exp() − 2 exp[ − (2J + D)]

m=1

+exp[ − ( + 4J + 2D)] − exp[( − D)] +exp[ − 2(J + D)]}={1 − exp[ − (2J + D + )]}2 :

(13)

We would like to remark that the above summation converges only for J + D +  ¡ 0. In two and three dimensions, the number of lattice animals for a given m grows exponentially, and we cannot 4nd a closed form for the partition function. We have only obtained results for m 6 15 for d = 2 and m 6 11 for d = 3, which was done by a straightforward application of the algorithm of Ref. [26]. Therefore, the evaluation of the total concentration Xt for the two- and three-dimensional cases presents some deviations, since it is a restricted sum in m of mXm . The deviation is very small in two dimensions, but can be appreciable in three dimensions. 4. Monte Carlo simulations We performed Monte Carlo simulation for the linear, square and cubic lattices in order to compare its results with those obtained by the independent-cluster approximation and the exact one-dimensional calculations. The simulations were carried out in the following way: (a) First, we established the model parameters (J , D, T , V and Xt ) and distributed, at random, VXt molecules in the lattice, each one occupying a single site. (b) We explored the phase space of the system by using the Kawasaki dynamics [29] (this dynamical process mimics the di5usion of particles in the system), where we choose, arbitrary, one amphiphile and try to exchange its position with one nearest-neighbor water (if there is not at least one water molecule, we restart the step (b)). The change is accepted or rejected following the Metropolis prescription [30]. One Monte Carlo step (1MCs) is de4ned as been V trials to exchange a pair of molecules. (c) After typically 105 MCs, the system reaches the thermal equilibrium and we begin to measure the quantities of interest. More 2 × 105 MCs are used to compute the mean values of the concentrations. The steps (a) – (c) are done for each sample and the error bars are given by standard deviation from 20 samples. We collected the data for di5erent temperatures and values of the ratio J=D in order to obtain the exponent
, the CMC curve and the phase diagram. The diagram in the plane J versus D exhibits the region where the system is in a micellized state. We must remember, however, that it is not a real phase diagram, since we have no true phase transition between the micellized and non-micellized states.

628

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

5. Results and conclusions In Fig. 2 we exhibit the CMC curves for the one-dimensional case. They present a plateau, which indicates that the amphiphiles added to the solution prefer to stay aggregated. The three full lines are exact results for di5erent temperatures and we also show the Monte Carlo results (circles) together with the independent cluster approximation (dotted curve). Note that the agreement between simulations and the exact curve is excellent, while the independent cluster approximation (ICA) overestimates the values of the concentration of the isolated molecules. This behaviour was expected because we have no cluster–cluster interaction, and the reaction c1 + cm  cm+1 (cm means a cluster of m molecules), that would reduce the concentration of isolated molecules, is not taken into account. Fig. 3 shows the aggregate-size distribution curve (d = 1 case) for two di5erent values of the parameter D, Xt = 0:1 (that is already above the CMC) and T = 0:45. The curve displays the local maximum and minimum, indicating we have a micellized system. For higher temperatures, this curve becomes monotonically decreasing and the system is not micellized. We call the di5erence between the vertical coordinates of the minimum and maximum concentrations as . We observed that, as D increases, bigger aggregates are favoured. Fig. 3 also exhibits the results of the simulations (circles) and the independent cluster approximation. Note that the ICA gives good quantitative results for the distribution of the aggregates. As a matter of

T=0.55 0.03 T=0.50

0.02

X1

T=0.45

0.01

0 0

0.05

0.1

0.15

Xt Fig. 2. Concentration of isolated amphiphiles X1 as a function of the total concentration Xt . The full lines are the exact results for the one-dimensional case and for three di5erent values of temperatures T (T = 0:45; 0:50 and 0.55). The circles are Monte Carlo simulations for T = 0:50 and the dotted line gives the independent cluster approximation at T = 0:50. Here J = 1 and D = 0:6. The dashed line is the ideal gas behaviour.

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

629

0.03

0.025

mX m

0.02

0.015 D=0.6

0.01

0.005 D=0.9

0 0

10

5

15

m Fig. 3. Volume fraction of aggregates mXm as a function of the aggregate size m for two values of the parameter D in d = 1. Full lines are the exact results for T = 0:45, Xt = 0:1 and J = 1. The dashed lines are the ICA calculations, and the circles are the results of simulations for D = 0:6.

2

1.5

βJ

M

1

0.5 NM

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

βD Fig. 4. Diagram in the plane J versus D obtained by exact calculations. Region M indicates the micellized states, while NM is the region where the system is not micellized. Here Xt = 0:1 in d = 1.

630

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633 0.006

0.005

∆

0.004

0.003

0.002

0.001

0 0.49

0.5

0.51

0.52

0.53

0.54

T Fig. 5. Di5erence in height  as a function of the temperature T . Exact results (continuous line) and simulations (circles) for Xt = 0:1, J = 1 and D = 0:6 in one dimension. The 4t of the simulation data to  ˙ (T − TM )
, gives
= 1:0 ± 0:1.

illustration, we show in Fig. 4 the region in the plane J versus D where the system is micellized for Xt = 0:1. It is not a true phase diagram because a phase transition between the micellized and non-micellized states does not occur. It is interesting to note that, as shown in Fig. 4, micellar aggregates appear only when both J and D are positive. It may be related to the fact that, in the one-dimensional case, the aggregation turns to be diKcult due to the spatial constraint, and both the terms in the hamiltonian must contribute to the aggregation. The di5erence in height between the vertical coordinates of the maximum and minimum in the aggregate-size distribution curve in one dimension is shown in Fig. 5 as a function of the temperature. The continuous line is the exact result while the circles represent the simulation data. For both data, the exponent is
∼ 1. In ICA the exponent is also
∼ 1, however, the micellization temperature TM is quite di5erent. We will now show the results for the two- and three-dimensional versions of the model. As we said before, only the Monte Carlo simulations and the ICA calculations are given. Fig. 6 shows the aggregation distribution curves for the model in two and three dimensions, and for a suitable choice of the parameters J and D. Again, the maximum and the minimum in these curves indicate that the system is micellized. In this 4gure we plotted the results of both approaches, simulations and ICA calculation, for a sake of comparison. As the dimensionality of the system increases, the ICA becomes less accurate, but their results are qualitatively sound. Note that for the twoand three-dimensional versions of the model, the parameters J and D must be positive and negative, respectively. This energetic competition provides the formation of

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

631

0.25

0.2

mX m

0.15

0.1

0.05

0 0

4

8

12

16

20

m

Fig. 6. Volume fraction mXm as a function of the aggregate size m for d = 2 and d = 3. Continuous line and connected triangles represent, respectively, the ICA results and the simulation data. Here, we have d = 2, T = 0:465, D = −0:7, J = 1 and Xt = 0:1. Dashed line and connected circles represent the ICA calculations and the simulation data, respectively. Results for d = 3, T = 0:49, D = −0:7, J = 1 and Xt = 0:1.

aggregates whose sizes are distributed according to the requirements to form micelles. This behavior is shown in Fig. 7, where we plot the diagram D versus J in two dimensions (the three-dimensional diagram is very similar). There we exhibit the region of the parameter space where micellar aggregates appear. Finally, Fig. 8 shows the plot of the di5erence in height  as a function of the temperature for d = 2 and 3. Again, the exponent
is near 1 independently of the spatial dimension. In conclusion, we have found exact results in one dimension, and performed Monte Carlo simulations in one, two and three dimensions for a lattice-gas model that mimics some characteristics of the micellar solutions. We have also applied an approximation scheme used for the imperfect gases (independent cluster approximation) to study the aggregation properties of the model. We observed that, for all spatial dimensions considered, the model captures the main characteristics of a micellized solution: the presence of a well-de4ned critical micellar concentration and the minimum and maximum in the aggregate-size distribution curve. The calculation of the exponent
, which measures the transition to a non-micellized state, furnishes the value 1, independently of the spatial dimension. We also determined the diagrams in the space of the parameters J and D where micellized states are possible. We veri4ed that the independent cluster approximation gives a good qualitative description of the model, but overestimates the concentration of isolated amphiphiles, indicating that the cluster–cluster interaction must play an important role in the process of aggregation as we increase the total concentration of the solution.

632

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633 0

-0.5

βD

M -1

-1.5

NM

-2 0.8

0.6

1

1.2

1.4

1.6

1.8

2

βJ Fig. 7. Diagram in the plane D versus J obtained by simulations. Region M indicates micellized states, while region NM is associated with non-micellized states. Here Xt = 0:1 in d = 2.

0.0025

(a)

0.0005

(b)

0.0004

0.0015

0.0003

0.001

0.0002

0.0005

0.0001

∆

0.002

0

0 0.5

0.502

0.504

0.528

0.53

0.532

T Fig. 8. Di5erence in height  as a function of the temperature T : (a) Simulation data (circles) for d = 2, D = −0:7, J = 1 and Xt = 0:1. (b) Simulation data (circles) for d = 3, D = −0:7, J = 1 and Xt = 0:1. The straight lines are the best 4t to the simulation data. For both dimensions we found
= 1:0 ± 0:1.

M. Girardi, W. Figueiredo / Physica A 324 (2003) 621 – 633

633

Acknowledgements This work was partially supported by the Brazilian agency CNPq. References [1] C. Tanford, Hydrophobic E5ect: 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] J. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1992. [4] A.Z. Panagiotopoulos, M.A. Floriano, S.K. Kumar, Langmuir 18 (2002) 2940. [5] U. Reimer, M. Wahab, P. Schiller, H.J. MSogel, Langmuir 17 (2001) 8444. [6] M. LITsal, C.K. Hall, K.E. Gubbins, A.Z. Panagiotopoulos, J. Chem. Phys. 116 (2002) 1171. [7] G. Cristobal, J. Rouch, J. CurIely, P. Panizza, Physica A 268 (1999) 50. [8] P.G. Bolhuis, D. Frenkel, Physica A 244 (1997) 45. [9] A.T. Bernardes, V.B. Henriques, P.M. Bisch, J. Chem. Phys. 101 (1994) 645. [10] J.N.B. de Moraes, W. Figueiredo, J. Chem. Phys. 110 (1999) 2264. [11] M. Girardi, W. Figueiredo, J. Chem. Phys. 112 (2000) 4833. [12] M. Girardi, W. Figueiredo, Phys. Rev. E 62 (2000) 8344. [13] M. Girardi, W. Figueiredo, Phys. Stat. Sol. A 187 (2001) 195. [14] M. Girardi, W. Figueiredo, Physica A 319 (2003) 421. [15] D. Duque, Phys. Rev. E 64 (2001) 063 601. [16] D. Duque, P. Tarazona, J. Chem. Phys. 107 (1997) 10 207. [17] T. Dalby, C.M. Care, Phys. Rev. E 59 (1999) 6152. [18] H. Wennerstrom, B. Lindman, Phys. Rep. 52 (1979) 1. [19] C.B.E. Guerin, I. Szleifer, Langmuir 15 (1999) 7901. [20] K.K. Chin, D.P. Landau, Phys. Rev. B 36 (1987) 275. [21] R.J. Baxter, F.Y. Wu, Phys. Rev. Lett. 31 (1973) 1294. [22] X.N. Wu, F.Y. Wu, J. Phys. A 22 (1989) L1031. [23] T.L. Hill, Statistical Mechanics, Principles and Selected Applications, Dover Publications, New York, 1987. [24] G. Gompper, M. Schick, Phys. Rev. B 41 (1990) 9148. [25] C.J. Thompson, Mathematical Statistical Mechanics, Princeton, NJ, 1972. [26] J. L. Martin, in: C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, Vol. 3, Academic Press, London, 1971. [27] S. Mertens, J. Stat. Phys. 58 (1989) 1095. [28] S. Mertens, M.E. Lautenbacher, J. Stat. Phys. 66 (1991) 669. [29] K. Kawasaki, in: C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, Vol. 4, Academic Press, London, 1974. [30] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087.