Obtaining pressure versus concentration phase diagrams in spin systems from Monte Carlo simulations
We propose an efficient procedure for determining phase diagrams of systems that are described by spin models. It consists of combining cluster algorithms with the method proposed by Sauerwein and de Oliveira where the grand canonical potential is obtained directly from the Monte Carlo simulation, without the necessity of performing numerical integrations. The cluster algorithm presented in this paper eliminates metastability in first order phase transitions allowing us to locate precisely the first-order transitions lines. We also produce a different technique for calculating the thermodynamic limit of quantities such as the magnetization whose infinite volume limit is not straightforward in first order phase transitions. As an application, we study the Andelman model for Langmuir monolayers made of chiral molecules that is equivalent to the Blume-Emery-Griffiths spin-1 model. We have obtained the phase diagrams in the case where the intermolecular forces favor interactions between enantiomers of the same type (homochiral interactions). In particular, we have determined diagrams in the surface pressure versus concentration plane which are more relevant from the experimental point of view and less usual in numerical studies.
pacs:05.10.Ln, 05.70.Fh, 05.50.+q
The importance of numerical simulations in Physics is due to the fact that very few models can be exactly solved. In principle one may directly simulate any model on a computer. Moreover, the Metropolis metr () and the Glauber glauber () algorithms used in Monte Carlo (MC) simulations are very general and easy to implement. In practice things are not so simple. Near second-order phase transitions the configurations generated by these algorithms present strong temporal correlations (critical slowing down), which prevent an efficient sampling of the configuration space. In addition, hysteresis effects due to metastability prevent a precise location of first-order transition lines.
In the last years, several techniques have been proposed to circumvent these problems such as the reweighting technique by Berg and Neuhaus berg () and the simulated tempering by Marinari and Parisi parisi (). A different approach is the use of cluster algorithms pioneered by Swendsen and Wang sw () and by Wolff wolff (). Several studies have shown the efficiency of cluster algorithms in reducing the critical slowing down sw (). More recently, it has been shown that the cluster dynamics may practically eliminate metastability in first order phase transitions So far this has been achieved only for the Blume-Emery-Griffiths (BEG) spin-1 model BEGMODEL (); BEG2 (), for which a special cluster algorithm has been developed bouabci (); rachadi ().
In this paper, we present a simple cluster algorithm that eliminates metastability in first-order phase transitions and combine it to the Sauerwein and de Oliveira (SO) method that allows us to obtain the surface pressure directly from numerical simulations sauerwein (). In the original formulation of the SO method, the authors used the Metropolis dynamics to generate the system configurations. However, as explained above, this is not the best choice near phase boundaries. We also introduce a simple procedure to calculate the order parameters and the concentrations of molecules in the neighborhood of a first-order line from numerical simulations.
As an application of our method, we have determined the phase diagrams of the Andelman model for Langmuir monolayers made of chiral molecules. More specifically, we were interested in surface pressure versus concentration phase diagrams, which are more interesting from the experimental point of view. The Andelman model was first studied by using the mean field approach on a bipartite lattice andelman1 (). However, X-ray diffraction experiments suggest that the condensed phases of Langmuir monolayers tend to form triangular structures, and not a bipartite lattice as considered by Andelman. In order to be more consistent with the physics of Langmuir monolayers, Pelizzola et al pelizzola () have studied the heterochiral case on a two dimensional triangular lattice, using the cluster variation method, and have obtained phase diagrams that are qualitatively different from Andelman’s. We study through MC simulations the remaining homochiral case, and we show that, in contrast with the heterochiral case, the MC and the mean field methods give results that are in good agreement.
This paper is organized as follows: in section 2 we briefly review the Andelman model, in section 3 we present the cluster algorithm and briefly review the Sauerwein de Oliveira method, in section 4 we discuss the numerical results and conclude in section 5.
Ii Chiral Langmuir monolayers and the Blume-Emery-Griffiths model
Langmuir monolayers are formed by spreading amphiphilic molecules in an air-water interface. Amphiphilic molecules are strongly asymmetric, constituted by two parts with opposite features. The first part—the head—, is hydrophilic. It is made of polar chemical groups and remains on the water. The second part—the tail—, is hydrophobic and made of hydrocarbon chains which remain in the air. When the tail is strongly hydrophobic, so that the molecules are insoluble in water, Langmuir monolayers form a quasi-two dimensional system. This system can be described in terms of the surface pressure and temperature and it displays several phases with different structural properties (see, for instance Ref kaganer ()).
A chiral molecule exists in two forms and , called enantiomers, related by a spatial transformation that involves a change of parity. An important feature of the physics of chiral Langmuir monolayers is the determination of the chiral discrimination, which occurs when the interaction energy between enantiomers with the same chirality is different from the interaction energy of enantiomers with different chirality. When the intermolecular forces favor the attraction between enantiomers of the same species, they are denominated homochiral and they lead to chiral segregation. On the contrary, if the attraction between different enantiomers is favored, they are named heterochiral and they lead to a racemic mixture.
To study the effect of the chirality theoretically, Andelman proposed a simple lattice gas model that can be described by the Hamiltonian pelizzola ()
where the first sum is over nearest-neighbor pairs, the letters and denote the sites of a two-dimensional triangular lattice, the letters and denote the enantiomer species ( or ), are the coupling energies ( and ), are the occupation numbers at site , and is the chemical potential of the species .
where . Thus, represents a () enantiomer and a vacancy. In this way, we obtain the BEG Hamiltonian
The case corresponds to the homochiral case (ferromagnetic BEG). When , we have the heterochiral one (antiferromagnetic BEG). We will concern ourselves with the homochiral case, since it has not been studied beyond the mean-field approach. The parameters depend on the interaction energies and between nearest-neighbor enantiomers through the formulae
The fields and are related to the chemical potential of the species and and they are given by
They are the conjugate parameters of the chiral order parameter and the density of enantiomers defined by
where the are the total number of enantiomers and is the number of lattice sites. In particular, we are interested in determining the concentration of enantiomers or ( or , respectively) given by
where , .
In this work we shall study first-order transitions between concentrated phases where the enantiomers are close to each other (phases and rich in enantiomers of type and , respectively) and the so called liquid expanded () phases, where there are many vacancies. In the spin-1 language, we shall study transitions between ferromagnetic and paramagnetic phases rich in zero spins.
Iii Monte Carlo method
iii.1 Cluster algorithm
In experiments involving chiral Langmuir monolayers the non chiral contribution to the interaction energy between the enantiomers is usually larger than the chiral one. In our simplified model, this corresponds to choosing the parameter larger than . In this paper, we will consider the ratio . This choice has also been previously made by Andelman andelman1 () and Pelizzola et al pelizzola (). For , Eq. (3) can be rewritten, up to a constant term, in the following way
where we used the following identity , , , , and is the coordination number. For this Hamiltonian, we propose the following cluster algorithm:
Choose randomly a site on the lattice and denote the value of its spin. This is the first spin of the cluster (seed).
Choose, with the probability 0.5, one of the two other possible spin values that are different from . Call this new value (it will remain fixed during the construction of the cluster). For example, if is , can be or .
Activate the links between the seed and its nearest neighbors that are equal to with probability . Each new spin connected to the cluster by an activated link is added to the cluster. Next, we repeat the activation procedure to all the new spins of the cluster. The process stops when all nearest neighbors have been tested and no new spin is accepted. Now, we attempt to change this cluster with spins equal to into a cluster with spins (see Fig. 1, for an example of a transition).
Evaluate the difference , where is the cluster bulk energy (calculated neglecting boundary links) when all spins are equal to and is the cluster bulk energy when all spins are equal to . If , we change all spins in the cluster to with probability . If , we change all spins in the cluster to with probability .
To prove that the algorithm satisfies the detailed balance condition we have to consider two types of transitions: and . For the first transition our algorithm is equivalent to Wolff’s wolff () and for this reason we shall concentrate on transitions of the second type. Let us consider, to exemplify, the transition . From Eq. (11), we obtain
where is the total number of boundary links that connect sites with spins inside the cluster and sites with spins outside the cluster.
The ratio between the transition probability and the reverse transition probability is given by
The bulk term is the sum of the probabilities associated with all possible ways of activating links to construct the cluster, the term is the probability of not including in the cluster a nearest neighbor site with occupation variable . Analogous comments hold for the transition . Clearly, because for each configuration of activated links in there is a corresponding one in (see Fig. 1). Recalling the definition of , given in step 4 of the algorithm, we see that the ratio of the flipping probabilities is always equal to . Finally, since the right hand sides of Eqs. (12) and (13) are equal and this equality implies detailed balance. It is worth mentioning that the algorithm proposed here is a particular case of the cluster algorithm proposed by Bouabci and Carneiro bouabci () and later extended by Rachadi and Benyoussef rachadi () for other regions of the parameter space.
iii.2 The Sauerwein and de Oliveira method
In order to determine the grand canonical potential from Monte Carlo simulations, one usually calculates one of its derivatives and numerically integrate the results. To use this technique one has to know the value of the grand canonical potential at a reference point and then numerically integrate along a path which connects the reference point to the point where one wants to calculate the grand canonical potential. An alternative is the method proposed by Sauerwein and de Oliveira sauerwein () that allows one to directly obtain the grand canonical potential from the MC simulation.
In this method, the largest eigenvalue of the transfer matrix is directly evaluated from Monte Carlo simulations. Since in the thermodynamic limit the grand partition function is proportional to the largest eigenvalue of the transfer matrix, its calculation enables us to determine all thermodynamic properties, in particular the surface pressure that is the negative of the grand canonical potential.
In order to explain how to obtain the largest eigenvalue, let us consider a triangular lattice with sites divided in successive layers with spins, (All this applies to the triangular lattice that we use in our paper.). The Hamiltonian may be decomposed in the following way
where due to the periodic boundary conditions . The probability of a given configuration of the system is given by
where is an element of the transfer matrix and
is the grand-canonical partition function. By using the spectral expansion of the matrix it is possible to show sauerwein () that
This expression enables us to calculate the largest eigenvalue of the transfer matrix in terms of the averages and , where when layers and are equal and zero otherwise. We use a MC simulation to generate the configurations with which we calculate the averages.
In the specific case of the BEG Hamiltonian in the triangular lattice with sites, the transfer matrix of a -layer is given by
The grand canonical potential per site in the lattice gas representation (or the free energy in the spin-1 representation) is given by
where is the surface pressure.
Iv Numerical results
In this section, we define the following dimensionless quantities:
where , the surface pressure, is given by Eq. (19).
As a check on the efficiency of the proposed cluster algorithm, we show in Fig. 2 the grand canonical potential versus the chemical potential for and . We considered a very low temperature, because in this case hysteresis effects are very strong. In Fig. 2 we compare the performances of the Metropolis and the cluster algorithm on a triangular lattice with periodic boundary conditions and linear dimension . To evaluate and to estimate its statistical error after equilibrating the systems we have used Monte Carlo steps divided into independent runs. Note that with the Metropolis algorithm the system is trapped in metastable states and even after millions of MC steps it does not undergo a transition to the stable phase. This does not happen with the cluster algorithm because the system is able to easily pass from one phase to the other. The efficiency of the algorithm allows us to determine first-order transition lines with high precision and the good quality of the data enables us to perform very precise finite size analysis.
In principle it is possible to determine the transition point using the free energy. As it can be seen in Fig. 2, there is a kink in the free energy as a function of at the transition point . One can then perform a finite size analysis to obtain . However, it is simpler and more efficient to analyze the susceptibility whose finite size behavior is well known for both first and second-order phase transitions. After determining we use the SO method at this point to calculate the surface pressure. In first-order phase transitions, the surface pressure, which is proportional to the negative of the grand canonical potential, does not have a finite size behavior as simple as the susceptibility. As a function of the system size, the surface pressure saturates quickly. Thus, the values of the surface pressure that we use in our graphs come from the largest lattices that we have simulated.
The susceptibility is defined as , where the magnetization . For a fixed system size , maintaining and fixed, and increasing towards the coexistence line, one observes a peak in the susceptibility at , as seen in Fig. 3, where the lines were drawn only to guide the eye. In the thermodynamic limit, this peak becomes a delta-function singularity. According to Refs. rBoKo (); challa (), the deviation of from its asymptotic value decays as , in agreement with our results, shown in the inset of Fig. 3 . From this law, we have obtained the extrapolated value .
In order to understand this result let us perform an exact zero temperature calculation of the transition point . At zero temperature the free energy , where is given in Eq. (3). The system chooses the phase that minimizes the energy . At , for small values of all spins are if or if . If is large enough, the energy is minimized when all spins are . The transition line is obtained by equating the energies in the ferromagnetic (condensed) phase with all (or ) and in the paramagnetic (liquid expanded) phase with all . Taking into account Eq. (3), with , and the definitions given in Eq. (20), we calculate the energy of the ferromagnetic phases and the energy of the paramagnetic phase
The equation gives the transition lines between the ferromagnetic and paramagnetic phases. The equation gives the transition line between the two ferromagnetic phases (this holds for , for the systems is in the paramagnetic phase). All these transition lines are represented in Fig. 4 that gives the phase diagram of the Langmuir monolayer in the plane of the chemical potentials (Recall that and .). The circles are the results of Monte Carlo simulations performed at . The error bars are smaller than the circles. It is interesting to remark that in the temperature interval relevant for Langmuir monolayers the zero temperature calculations give practically the same results as the MC simulations and the mean-field calculations for the transition lines. In Langmuir monolayers language, for and low values of (higher chemical potentials), we have the condensed phase characterized by a mixture of the two enantiomers. For higher values of a transition from the condensed phase, rich in enantiomers , to the phase poor in enantiomers, the liquid expanded phase, takes place. When the chemical potential of the species are different (), we have larger fraction of enantiomers whenever , and in the limit of the solution only contains the enantiomer .
Another procedure for locating the phase transition consists in determining the crossing point of the versus isotherms for different system sizes. As showed in the Ref. fiore4 (), the crossing point is independent of the lattice size and properly identifies the transition, as shown in Fig. 5. We shall present below an independent derivation of this important result based on the work of Borgs and Kotecký rBoKo ().
For all curves of versus cross at and for this value of . This criterion for estimating the value of for which the phase transition takes place agrees very well with the finite size analysis of the susceptibility that we have discussed above. For , two phases coexist at the point which now depends on and all isotherms cross at . We remark that if single flip algorithms are used to generate the dynamics, one will not be able to determine the crossing of the curves due to hysteresis effects.
More relevant from the point of view of Langmuir monolayers, and other physical systems involving mixture of molecules, are the surface-pressure versus concentration diagrams. But before discussing this phase diagram we will describe our procedure to fit the curves in Fig. 5 and to obtain the limit of and that are used to determine the concentration (see Eq. (10)). Since the simulated system is finite, the calculated quantities will be affected by finite size effects. As mentioned previously, in the last years, the finite size theory of first order phase transitions has been studied extensively for quantities, such as the specific heat and the susceptibility. There are fewer studies for the dependence on the system size of quantities like the magnetization or the concentration of molecules gupta (); borgs94 (). In the following, we propose a method to determine the concentrations of the phases that coexist directly from the numerical simulations. The first step consists in noting that (or ) isotherms can be fitted by the equation
where , , and are fitting parameters and . We are going to show below that depends on the system size and the temperature . An analogous expression can be written down for the order parameter . The expression above was inspired by the work of Borgs and Kotecký rBoKo (), where it is shown that at low temperatures the partition function for two coexisting phases can be written as
where is the magnetic field (our system also depends on the crystal field ), is a constant of the order of the infinite volume correlation length and is a metastable free energy for the phase ().
In our system, for three phases coexist at the triple point (). Thus, we expect the sum of three exponentials instead of two as in Eq. (23). We have assumed that all three exponentials have the same weight (we shall use our results to check this point). In the neighborhood of the triple point
where the , , are respectively the metastable free energies of the paramagnetic and ferromagnet phases. Away from the coexistence curve, only the associated with the correct phase remains and becomes the free energy of the system ( ).
The parameters and are given by
At the triple point , where , and the exponentials in Eq. (26), which contain the only dependence on the lattice size, cancel out and we obtain
This is the reason why the curves for different lattice sizes cross at the same point. The crossing point provides another method do locate the phase boundaries.
Our calculations are performed at low temperatures. In mean-field, the temperatures are measured in units of the coordination number . For the triangular lattice, . Our MC temperature is equivalent to a mean-field temperature. We can use the exact zero temperature energies given in Eq. (21) to estimate the derivatives in Eq. (27). Recalling that , , at , we obtain and . Thus, at the triple point which is the result that we obtain in Fig. 5.
In the curves plotted in Fig. 5, varies in the interval which is very narrow. It is possible in this case to expand the , , around the triple point,
where , and , for ,
which has the same form as Eq. (22) after we divide the numerator and the denominator by .
In Fig. 5 the symbols stand for the values of obtained from the numerical simulations and the solid lines are fits of the points using Eq. (22) by minimizing the merit function press (). In order to perform the fittings we used the Levenberg-Marquardt method that is well described in Ref. press (), where one can also find the subroutines that are necessary to implement the method. These subroutines return the variances of the fitting parameters and the quality of the fitting. A few words about the implementation of the subroutines is in order. Our fitting function Eq. (22) contains exponentials whose arguments may become very large. In order to avoid numerical overflow it is convenient to use the asymptotic values of when becomes too large. Define, for example, for and for . Of course, the number is rather arbitrary. Non-linear fittings depend on a good initial guess of the fitting parameters. One may proceed as follows. Note that and . Since in the simulations the interval is finite, instead of taking the limit we use the values of in our data set associated with the largest and the smallest values of . Call them and , respectively and put , . Next define and , where , is chosen in the region where the graph has already started to curve down. It is simple to solve the fitting parameters in terms of these quantities.
With this choice for the initial parameters, the convergence of the fitting routine is very fast and the quality of the fitting is very good (the factor that measures the goodness-of-fit press () is close to ). In Tables 1 and 2 the errors of the parameters are the square roots of the variances (standard deviations) that are returned by the fitting routines.
The fitting parameters for the curves in Fig. 5 are given in Table 1. Now we can check the equal weight hypothesis for the exponentials. For the two condensed phases have the same free energy and two of the three exponentials are identical, as we discussed above. The curves in Fig. 5 are in the vicinity of the triple point, so we expect that . This is the result that we obtain (see Tables 1 and 2 for the magnetization ).
For , there is the coexistence of two phases ( and or ). Near the transition we have the sum of two exponentials with the same weight, as in Eq. (23). We have checked that near the transition line, as it was expected.
Comparing Eqs. (22) and (26) we note that the parameter that appears in the exponent should be proportional to the system volume. A plot gives the straight line with and for table 1; and and for table 2. Thus, as expected, the constant scales with the volume of the system.
The values of for the condensed and liquid expanded phases are calculated by taking the limit in Eq. (22). The condensed phases occur in the region and for this reason as . The liquid expanded phase occurs in the region and as . The curves for are very similar to the curves in Fig. 5 and can be also be fitted by an expression analogous to Eq. 22. Having calculated and , we obtain through the expression . The use of instead of is due to technical reasons (see sections 2.3.3 and 2.3.4 in Ref. binder ()). As a consequence, the magnetization is small but not zero when . This introduces a small distortion in the diagram of Fig. 6 near , but symmetry arguments guarantee that when and the coexistence curve passes through the point with (filled circle in Fig. 6).
The half side of the diagram is mapped onto the right hand side of the surface-pressure versus concentration diagram () whereas corresponds to the concentration range. As mentioned above, implies that the fraction of enantiomers and are equal and in the coexistence of the three phases, one has . From the point of view of homochiral Langmuir monolayers, the chiral segregation takes place, in contrast to the heterochiral case, in which one has a racemic mixture. The surface pressure of a 1:1 mixture of enantiomers is higher than the pressure for pure enantiomers. This feature is verified in experiments in which the chiral segregation occurs kaganer (). Unfortunately, to date few experiments have been performed covering the whole range of concentrations, usually they are restricted to the 1:1 mixture and the pure cases. For comparison, we have also plotted in Fig. 6 the results obtained from the mean field technique.
In contrast to the heterochiral case, for which the mean field results disagree with those obtained from the cluster variational method, in the homochiral case the accordance between mean field and the numerical simulations is very good.
In this paper, we present an efficient way for determining phase diagrams from numerical simulations. To illustrate it, we have considered a simple model that describes the behavior of homochiral Langmuir monolayers, which is equivalent to the BEG model. It is worth mentioning that although we have interpreted the phase diagrams obtained here in terms of Langmuir monolayers, similar phase diagrams are obtained when one uses the BEG model to describe a mixture of two distinct species with vacancies. The use of a cluster algorithm that eliminates metastability in first order phase transitions allows us to precisely locate the first-order transitions lines. To determine the surface pressure we used the method proposed by Sauerwein and de Oliveira in which the surface pressure is determined directly from the numerical simulations without the necessity of performing numerical integrations. The fitting procedure proposed in this paper to determine the concentrations, based on the work of Borgs and Kotecký rBoKo (), is easy to implement and uses all information contained in the order parameter curve. It seems to improve on the usual finite size analysis for the magnetization near first-order transition lines in that it does not present “overshooting” effects gupta (); borgs94 () and both and present a monotonic behavior as a function of , but this point has to be further investigated by increasing the statistics. The elimination of metastability also enables us to use the crossing of the curves (or ) for different lattice sizes as a criterium for locating the phase boundaries. This usually cannot be done due to hysteresis effects. Finally, we remark that our approach is general and it can be used for any spin model. In systems for which a cluster algorithm is not available, we can use other techniques, such as the multicanonical approach berg () or the simulated tempering parisi () to generate the dynamics.
C. E. F. acknowledges the financial support from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Grant No. 06/51286-8.
- (1) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
- (2) R. J. Glauber, J. Math. Phys. 4, 294 (1963).
- (3) B. A. Berg and T. Neuhaus, Phys. Lett. B 267, 249 (1991); Phys. Rev. Lett. 68, 9 (1992).
- (4) E. Marinari and G. Parisi, Europhys. Lett. 19(6), 451 (1992).
- (5) R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58, 86 (1987).
- (6) U. Wolff, Phys. Rev. Lett 62, 361 (1989).
- (7) M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev. A 4 , 1071 (1971).
- (8) W. Hoston and A. N. Berker, Phys. Rev. Lett. 67, 1027 (1991).
- (9) M. B. Bouabci and C. E. I. Carneiro, Phys. Rev. B 54, 359 (1996).
- (10) A. Rachadi and A. Benyoussef, Phys. Rev. B 68, 064113 (2003).
- (11) R. A. Sauerwein and M. J. de Oliveira, Phys. Rev. B, 52, 3060 (1995).
- (12) D. Andelman, J. Am. Chem. Soc, 111, 6536 (1989).
- (13) A. Pelizzola, M. Pretti and E. Escalas, J. of Chem. Phys. 112, 8126 (2000).
- (14) C. Borgs and R. Kotecký, J. Stat. Phys. 61 79 (1990) and Phys. Rev. Lett. 68 1734 (1992).
- (15) M. S. S. Challa, D. P. Landau and K. Binder, Phys. Rev. B 34, 1841 (1986).
- (16) V. M. Kaganer, H. Möhwald, and P. Dutta, Rev. Mod. Phys. 71, 779 (1999).
- (17) C. E. Fiore, V. B. Henriques and M. J. de Oliveira, J. Chem. Phys. 125, 164509 (2006).
- (18) S. Gupta, A. Irbäck and M. Ohlsson, N. Phys. 409, 663 (1993), S. Gupta and A. Irbäck, N. Phys. 30, 861 (1993).
- (19) C. Borgs, P. E. L. Rakow and S. Kappler, J. of Phys. I, 4, 1027 (1994).
- (20) W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge Univ. Press, New York, 1987).
- (21) K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics (Springer-Verlag, New York Berlin Heidelberg, 1992).