Investigation of Fluid–Fluid and Solid–Solid Phase Separation of

Article pubs.acs.org/Langmuir

Investigation of Fluid−Fluid and Solid−Solid Phase Separation of Symmetric Nonadditive Hard Spheres at High Density W. T. Gózd́ ź* Institute of Physical Chemistry Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland ABSTRACT: We have calculated the values of the critical packing fractions for the mixtures of symmetric nonadditive hard spheres at high densities for small values of the nonadditivity parameter. Calculations have been performed for solid−solid and fluid−fluid demixing transitions. A cluster algorithm for Monte Carlo simulations in a semigrand ensemble was used, and the waste recycling method was applied to improve the accuracy of the calculations. The finite size scaling analysis was employed to compute the critical packing fractions for infinite systems with high accuracy.

■

INTRODUCTION Symmetric binary mixtures are a very convenient model to study various physical phenomena. In such mixtures, it is assumed that the interaction between components A is the same as the interaction between components B and only the interaction between components A and B is different. In the Lennard−Jones model, this difference may be introduced by setting different values of the depth of the potential well ϵ for the interaction of different species. For example, the choice ϵAB = 0.65ϵAB will lead to liquid liquid separation. Such a model was studied with great success by Keith Gubbins in the context of liquid mixtures in porous materials.1−4 Symmetric mixtures of nonadditive hard spheres behave in a similar manner. The behavior of the hard core fluid mixtures has been studied in bulk5−12 and in restricted geometry13−16 in three-dimensional and in two-dimensional systems.11,16−21 So far the investigation of the symmetric mixtures was mainly focused on the fluid− fluid phase transition. The simulations of the solid−solid transition for the nonadditive hard spheres have not been performed. Here we would like to concentrate on the behavior of the mixtures at high densities since they have not been yet examined at such conditions. Moreover the proximity of the fluid−solid phase transition makes the behavior of the system very interesting. We study the mixture of symmetric nonadditive hard spheres with the interaction potential defined by ⎧ r < σαγ ⎪∞ if Uαγ(r ) = ⎨ ⎪ ⎩ 0 if r > σαγ

where Δ is the nonadditivity parameter. We study the mixtures of positive nonadditivity with different values of the parameter Δ. This is evidently an athermal model, because all allowed configurations are of the same energy. Such mixtures can nevertheless separate into two phases: one phase rich in the component A and the other one rich in the component B. The phase separation is induced by the competition between the entropy of mixing and the entropy associated with the volume available for the particles. The increase of the packing fraction defined as η = ηA + ηB =

© XXXX American Chemical Society

(3)

where V is the volume of the system and NA and NB are the numbers of A and B particles, leads to larger decrease of the available volume in the homogeneous mixture than in the phase-separated system. This effect starts to dominate over the decrease of the entropy of mixing in the phase-separated system for η > ηc . Thus, ηc plays a role analogous to the critical temperature Tc , and (ηc − η)/ηc plays a role analogous to (T − Tc)/Tc in non-athermal mixtures. Simulation results10,11 show that this model system belongs to the Ising universality class. At high densities in a one component hard spheres system, a fluid−solid phase transition takes place. We may expect that the fluid−solid phase transition will be also present for the mixtures of nonadditive hard spheres if the nonadditivity parameter is sufficiently small. It may be also expected that for a certain range of the density there will be a region of coexistence of

(1)

where r is the distance between the centers of two spheres, indexes α ∈ {A , B} and γ ∈ {A , B} describe the species. The length scale is set by the A-component hard-sphere diameter, σAA = 1. For symmetric nonadditive mixtures, σBB = σAA and σAB = (σAA + σBB)(1 + Δ)/2

3 3 π NAσAA π NBσBB + 6 V 6 V

Special Issue: Tribute to Keith Gubbins, Pioneer in the Theory of Liquids Received: June 28, 2017 Revised: August 7, 2017 Published: August 21, 2017

(2) A

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir solid and fluid phases. Here, I concentrate on the regions of the density where there are only two coexisting solid or fluid phases. In the absence of an external field, the symmetry of interactions implies for the two coexisting phases I and II the following relations: xAI = x BII ,

x BI = xAII

(4)

μAI = μAII = μBI = μBII

(5)

and

where μαi and xαi are the chemical potential and the composition of the component α in the ith phase, respectively.1,10 Here the difference of the chemical potentials h = μA − μB plays the role of an external field. We are interested in the critical properties of the mixture when the external field is zero. Along the symmetry line h = 0, the composition of the coexisting phases is symmetric; thus, the critical point is at the concentration xc = NA /(NA + NB) = 0.5. We model an open system in contact with a reservoir by using the semigrand10,22 Monte Carlo23,24 simulation method. In this method, the system is simulated under a constant total number of particles N, total volume V, temperature T, and the difference of the chemical potential of one species with respect to an arbitrarily chosen species Δμ. Thus, the number of molecules of each species fluctuates, while the total density remains constant. The semigrand ensemble is superior to the grand ensemble in simulating dense fluids, because there are no particle insertion moves which are inefficient for dense fluids. For a symmetric binary hard spheres mixture, the internal energy and the chemical potential difference between the components of the mixture are both zero. The Monte Carlo simulations in the semigrand ensemble Monte Carlo simulation consist of two kinds of moves: translation and identity change. The identity change moves can be implemented in the following way. A molecule can be chosen randomly from all the molecules and the identity change move is always accepted if there is no overlap between the particles after the identity change. The calculations of the solid−solid phase equilibria in standard grand canonical ensemble or even in semigrand ensemble are inefficient due to low acceptance ratio of the Monte Carlo moves. Therefore, we use a cluster algorithm to perform the simulations.7,16,21 The idea of the cluster moves is the following. The system is divided into a set of clusters. The molecules belong to a cluster if the distance from any molecule to any other molecule is smaller than σAA + Δ. The schematic representations of the clusters build with such a procedure is presented in Figure 1. When the clusters are formed, the identity of all molecules in each cluster is randomly changed with the probability p = 0.5. We identify the clusters using the SLINK algorithm,25−27 which is fast and does not require huge amount of computer memory. In the first step of the single-linkage clustering, each molecule forms a cluster of its own. At each step, two clusters separated by the shortest distance are combined. The process ends when all the molecules form one large cluster. At the end we obtain the dendrogram which shows which clusters are combined at a given step. From that dendrogram we construct a set of clusters for a given maximal separation, as shown in Figure 1. We have performed the calculations using only cluster moves because for

Figure 1. Schematic representation of the clusters of particles obtained in SLINK procedure. The blue an red circles represent the particles of different type. The yellow ring has the width equal to Δ/2 .

the dense systems single identity change moves become inefficient due to a low acceptance ratio. For the translation moves, the maximum displacement is chosen to obtain 30% acceptance ratio. The identity-change and the translation moves are chosen randomly with the ratio of N translations per one cluster identity change move. We have performed the calculations in a cubic simulation box with periodic boundary conditions. The initial configuration was set up on a face centered (fcc) or simple cubic (sc) lattice. For the solid−solid phase separation, we used fcc latice, for the fluid− fluid phase separation we used mainly sc lattice. The averages are taken over 107 Monte Carlo cycles, where the cycle consists of N translations and one cluster identity change move. During the simulation, we obtain the probability density function of the concentration of the component A for a given number of molecules PN (x), where x = NA /(NA + NB). We use the waste recycling method28 to reduce the errors in the calculations of PN (x). The procedure which implements the waste recycling method is described below. Let us define the order parameter M as M = NA − NB

(6)

Let us denote the difference in the number of particles of each kind in a cluster j by mj The difference between all particles of each kind in the system with the total number of particles N = NA + NB is thus n

M=

∑ mj j=1

(7)

We can compute the characteristic function for such defined order parameter as f (k) = exp(ikM )

(8)

where the angular brackets denote the thermal average and i denotes the imaginary number i 2 = −1. The characteristic function is defined as the Fourier transform of the probability density function. Here, it is the probability density of the order parameter M. Equation 7 can be used to transform the exponential function in eq 8 in the following way B

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

in a hard spheres system from fluid to solid.29 We can expect that for the mixture of nonadditive hard spheres similar fluid to solid transition is present. For large nonadditivity parameter and large density, the mixture will separate into two phases pure or almost pure phases and these phases will behave as a one component hard sphere system. For small nonadditivity parameters, the repulsion between unlike molecules takes place at such a short distance that the behavior of the mixture is very similar to the behavior of the pure fluid with respect to fluid− solid phase transition. There will be a coexistence region between the fluid and solid phases for some range of the packing fraction. The number of fluid and solid phases at coexistence depends on the value of the nonadditivity parameter Δ. If Δ is small, the mixture separates only in the solid phase and the fluid and solid phases at coexistence will have the concentration xA = 0.5. If Δ is large, then the mixture separates already in the fluid phase and is separated also in the solid state. In this case, two fluid phases will coexist with two solid phases. We have observed that at some range of Δ the mixture which is separated at the fluid phase is mixed again at low densities of the solid phase and separates for higher densities of the solid phase. The typical configurations for the solid and fluid phase for the hard spheres mixture are shown in Figure 3. The solid phase is the fcc crystal. We have also

n

exp(ikM ) =

∏ exp(ikmj) j=1

(9)

If the system is divided into n clusters, we can construct 2n submoves since each cluster can have two values of the order parameter either mj or −mj . We can take all these configurations into account and calculate the submove average for the total cluster move c. This submove average can be denoted by fc and calculated in the following way. n

fc =

∏ cos(kmj) (10)

j=1

Without the waste recycling, only one configuration is taken to calculate averages. When the waste recycling method is used at each cluster move, 2n configurations are taken to calculate averages, where n is the number of clusters. Therefore, the statistical errors are much smaller. f (k) can be estimated by the following equation C

fest (k) =

1 ∑f C c=1 c

(11)

where C is the total number of cluster moves. The probability density P(M) to find the system with a given value of the order parameter M can be obtained by Fourier transform of eq 11. The P(M) can be easily transformed into different probability distribution functions more frequently used in the description of binary mixtures. If we replace M by NA = (M + N )/2 and divide NA by the total number of particles N = NA + NB, we obtain the probability density of concentration P(x). This is just a change of variable from M to x. Figure 2 shows the probability density of the concentration x = NA /(NA + NB) calculated without and with the help of the waste-recycling method.

Figure 3. Left side: fluid phase η = 0.4687 , N = 63 × 4 , and Δ = 0.04 . Right side: fcc solid phase η = 0.5300, N = 63 × 4 , and Δ = 0.04 .

performed the calculations for simple cubic (sc) lattice and body centered cubic (bcc) lattices. The sc lattice is unstable and the sc crystal never forms. The bcc crystal is formed but in the crystal phase the mixture is always separated into two almost pure one component phases. This can be understood if we compare the closest packing for these lattices which is 0.74 for fcc and 0.68 for bcc. In molecular simulations the results of calculations depend on a system size. The dependence is more pronounced for the calculations near the critical point since the correlation length becomes larger and larger when the system is closer and closer to the critical point. When the correlation length is larger than the size of the computational box, the results of the calculations become biased. This is why the exact calculations of the critical point parameters may be very difficult to perform. When the system is in the two phase region, we do not always have only one phase in the simulation box during the simulation in the semigrand Monte Carlo method. It is possible that in the simulation box we will have either the first or the second phase. With some frequency, the first phase disappears and the second phase appears and vice versa. The frequency of this change is higher the closer the system is to the critical

Figure 2. Comparison of the histograms of the concentration PN (x) calculated with the help of the waste recycling method (red curve) and just with the cluster algorithm (black curve) for Δ = 0.04 , η = 0.4686, and N = 63 × 4 .

■

RESULTS AND DISCUSSION We have calculated the critical packing fractions ηc for the range of small values of the nonadditivity parameter Δ. In the former studies, the calculations have been performed for Δ ≥ 0.1, where only fluid−fluid phase separation is present. For the calculations for Δ ≤ 0.1, we have a problem with efficiency of the simulation method due to high density. Moreover, at the packing fraction close to 0.5, there is a phase transition present C

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir point. Because of that, it is hard to calculate the concentration at the coexistence as an ensemble average and to obtain the location of the critical point from the shape of the coexistence curve. One may try to overcome this problem by taking the most probable value of the concentration instead of the ensemble average, but near the critical point this procedure may be problematic. The critical exponent for the three-dimensional fluids is β ≈ 0.3258. The distribution of the concentration is not sharply peaked and it is difficult to obtain the most probable value of the concentration with a sufficiently high accuracy. Fortunately, it is possible to use the calculations performed for the systems of finite size to extract the information about the infinite systems by using the concept of the finite-size scaling.10,30−32 For each value of the nonadditivity parameter Δ, the critical packing fraction of the infinite system, ηc (∞), can be calculated from the set of apparent critical packing fractions in systems with N particles, ηc*(N ), according to the relation:10,31,32 ηc*(N ) − ηc (∞) ∝ N −(1 + θ)/(dν)

spheres, compared with the distribution function for the 3D Ising model. PN*(y) agrees very well with the universal distribution P*(y) for the value of η that we identify with the apparent critical packing fraction ηc*(N ). In practice, the shape of the rescaled distribution function PN*(y) is compared with the universal distribution of the order parameter, P*(y), for the first estimation of the apparent critical packing fraction. In order to determine the precise value of the apparent critical packing fraction, the fourth order cumulants, UN = 1 −

M4 3 M2

2

(13)

are calculated for the values of the packing fraction η for which the distribution functions are the most similar to the universal distribution P*(y). To calculate the matching point, the values of the cumulant UN (η) have been interpolated near the universal value U N* . We use the fourth order cumulant UN ,33 because its value at the fixed point has been already calculated for the 3D Ising universality class. The nth moment M n can be easily calculated from the distribution of the order parameter PN (M ):

(12)

where the critical exponent ν = 0.629 for the 3D Ising universality class and d = 3 is the space dimension. The apparent critical packing fraction for a finite-size system can be determined from the distribution of the order parameter, PN (M ), calculated in the simulations as a histogram. M is the order parameter defined as M = NA − NB. The distribution PN (M ) is rescaled in such a way to have unit norm and unit variance, and the rescaled distribution is denoted by PN*(y), where the rescaled order parameter is y = aM −1N β /(dν)M with aM denoting a proportionality constant. For the apparent critical packing fraction ηc*(N ), PN*(y) has a universal shape P*(y).10,32 Since this model belongs to the Ising universality class, the universal function P*(y) is known.32 The rescaled function PN*(y) is compared with P*(y) for several values of η, and the best fit gives us ηc*(N ). Figure 4 shows the orderparameter distribution function PN*(y) calculated in the simulations for the system with Δ = 0.1 and N = 43 × 4

n

⟨M ⟩ =

∑M M nPN (M ) ∑M PN (M )

(14)

The sum is taken over all values of the order parameter M, where M = NA − NB. The apparent critical packing fractions ηc*(N ) satisfy the equation UN (ηc*) = U *, and have been read off from the interpolated line for the value of U * ≈ 0.47 .34 In Figure 5, we present the extrapolation of the apparent critical packing fractions ηc*(N ) calculated for finite values of N to the critical packing fraction for the infinite system ηc . The values of the critical packing fraction ηc for the nonadditivity parameter Δ ≤ 0.1 are plotted in Figure 6 and listed in the Table 1. The value of ηc for Δ = 0.1 is in good agreement with previous calculations.10 In Figure 6, we see two

Figure 5. Extrapolation of the critical packing fraction to the thermodynamic limit from the apparent critical packing fractions estimated by matching the order-parameter profile PN (M ) for a few system sizes N = 6 3 , 7 3 , 8 3 , N = 53 × 2, 63 × 2, 73 × 2 , N = 43 × 4, 53 × 4, 63 × 4 . The infinite-volume critical packing fraction from extrapolation to N = ∞ is ηc = 0.330736(20). θ = 0.54 . ν = 0.629. Δ = 0.1.

Figure 4. Normalized distribution of the order parameter PN (y) (open circles) for N = 43 × 4 system size, expressed as the function of the scaling variable y = aM −1N β /(3ν)M , at the packing fraction η = 0.3295, Δ = 0.1. The full curve is the fixed-point function for the Ising model.32 ν = 0.629. D

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

out that at some small range of the nonadditivity parameter this increase of the entropy of solid is so large that the mixture can again become totally mixed despite the increase of the density and related to that the decrease of the free volume available to the molecules. For small values of Δ, there is no phase separation in the fluid and such phenomenon cannot happen. For large values of Δ in the solid phase, the decrease of the entropy related to the free volume would be too high to allow for total mixing of the components, and therefore, it is preferable to separate into two solid phases with the concentration different than xA = 0.5.

■

SUMMARY We have determined for the first time the critical packing fractions for the symmetric nonadditive hard spheres model for low values of the nonadditivity parameter Δ ≤ 0.1. The phase transitions for the solid phases were examined for the first time. We have discovered the phenomenon of mixing the separated mixture after the phase transition to the solid phase. We have shown how the waste recycling method can be used together with the cluster algorithm to increase the accuracy of calculations of probability density functions.

Figure 6. Critical packing fraction ηc as a function of the nonadditivity parameter Δ. The dashed lines denote the values of the packing fraction determined for the hard spheres fluid and solid phase at coexistence. The solid squares denote the value of the packing fraction for the closest packed fcc lattice.

Table 1. Critical Packing Fractions ηc for Infinite Systems for Different Values of the Nonadditivity Parameter Δ for Fluid and Solid Phases fluid

■

solid

Δ

ηc

Δ

ηc

0.1 0.09 0.08 0.07 0.06 0.05 0.045 0.04 0.035

0.330736(20) 0.347050(20) 0.365244(35) 0.385842(90) 0.409333(150) 0.436751(250) 0.452416(300) 0.469639(500) 0.488552(550)

0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005

0.532210(200) 0.555441(50) 0.579056(50) 0.603295(100) 0.628417(100) 0.654585(100) 0.681877(50) 0.710483(100)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

W. T. Gózd́ ź: 0000-0003-4506-6831 Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS Support from NCN Grant No. 2015/19/B/ST3/03122 is acknowledged.

branches of points, one for solid−solid (open squares) and one for liquid−liquid (open circles) phase separation. For low values of Δ, the phase separation is only in the solid state. This is the separation into two solid phase one rich in the component A and the other rich in the component B. When Δ is larger, we naturally observe the fluid−fluid phase separation. The most interesting however is the behavior of the mixture for a range of the nonadditivity parameter where we see two critical points at different densities for the same Δ. This would indicate that the mixture which is separated in the fluid phases is mixed when the solid phase is formed at the higher packing fraction. When the packing fraction is increased, the mixture separates again into two solid phases. Since we are dealing with an athermal system with zero energy, all these phase transition must be governed by the entropy. This is either the entropy related to the free volume available to the molecules or the entropy of mixing. At low density, the increase in the entropy of mixing is bigger than the decrease of the entropy related to the free volume. This leads to the complete mixing of the two component hard spheres fluid. When the density is increased, this is no longer the case and the mixture separates into two phases. One may expect that the increase of the density for fixed Δ will always lead to the demixing of the nonadditive hard spheres mixtures. However, as it is shown in Figure 6, this is not always the case. In a one component system, the transition from fluid to solid phase takes place because the entropy of the solid is larger than the entropy of the fluid at the same density. It turns

■

REFERENCES

(1) Gózd́ ź, W. T.; Gubbins, K. E.; Panagiotopoulos, A. Z. Liquidliquid phase transitions in pores. Mol. Phys. 1995, 84, 825−834. (2) Gelb, L. D.; Gubbins, K. E. Liquid-liquid phase separation in cylindrical pores: Quench molecular dynamics and Monte Carlo simulations. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1997, 56, 3185−3196. (3) Gelb, L. D.; Gubbins, K. E. Kinetics of liquid-liquid phase separation of a binary mixture in cylindrical pores. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1997, 55, R1290−R1293. (4) Gelb, L. D.; Gubbins, K. Studies of binary liquid mixtures in cylindrical pores: phase separation, wetting and finite-size effects from Monte Carlo simulations. Phys. A 1997, 244, 112−123. (5) Amar, J. G. Application of the Gibbs ensemble to the study of fluid-fluid phase equilibrium in a binary mixture of symmetric nonadditive hard spheres. Mol. Phys. 1989, 67, 739−745. (6) Lomba, E.; Alvarez, M.; Lee, L.; Almarza, N. Phase stability of binary non-additive hard-sphere mixtures: A self-consistent integral equation study. J. Chem. Phys. 1996, 104, 4180−4188. (7) Johnson, G.; Gould, H.; Machta, J.; Chayes, L. Monte Carlo study of the Widom-Rowlinson fluid using cluster methods. Phys. Rev. Lett. 1997, 79, 2612−2615. (8) Saija, F.; Pastore, G.; Giaquinta, P. V. Entropy and Fluid-Fluid Separation in Nonadditive Hard-Sphere Mixture. J. Phys. Chem. B 1998, 102, 10368−10371. (9) Gózd́ ź, W. T. Phase separation of binary mixtures of symmetric non-additive hard spheres. Polym. J. Chem. 2001, 75, 517−525. E

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir (10) Gózd́ ź, W. T. Critical-point and coexistence curve properties of a symmetric mixture of nonadditive hard spheres: A finite size scaling study. J. Chem. Phys. 2003, 119, 3309−3315. (11) Buhot, A. Cluster algorithm for nonadditive hard-core mixtures. J. Chem. Phys. 2005, 122, 024105. (12) Pellicane, G.; Pandaram, O. D. Gibbs ensemble Monte Carlo of nonadditive hard-sphere mixtures. J. Chem. Phys. 2014, 141, 044508. (13) Duda, Y.; Pizio, O.; Sokolowski, S. Nonadditive Binary Hard Sphere Mixture in Disordered Hard Sphere Matrices: Integral Equations and Computer Simulation. J. Phys. Chem. B 2004, 108, 19442−19450. (14) Gózd́ ź, W. T. Phase separation of the Widom−Rowlinson mixture in confined geometry. J. Chem. Phys. 2005, 122, 074505. (15) De Sanctis Lucentini, P. G.; Pellicane, G. Critical Behavior of Symmetrical Fluid Mixtures in Random Pores. Phys. Rev. Lett. 2008, 101, 246101. (16) Almarza, N. G.; Martín, C.; Lomba, E.; Bores, C. Demixing and confinement of non-additive hard-sphere mixtures in slit pores. J. Chem. Phys. 2015, 142, 014702. (17) Fiumara, G.; Pandaram, O. D.; Pellicane, G.; Saija, F. Theoretical and computer simulation study of phase coexistence of nonadditive hard-disk mixtures. J. Chem. Phys. 2014, 141, 214508. (18) Muñ oz-Salazar, L.; Odriozola, G. Phase behaviour and separation kinetics of symmetric non-additive hard discs. Mol. Simul. 2010, 36, 175−185. (19) Nielaba, P. Phase transitions and quantum effects in adsorbed monolayers. Int. J. Thermophys. 1996, 17, 157−167. (20) Saija, F.; Giaquinta, P. V. Monte Carlo simulation and phase behavior of nonadditive hard-core mixtures in two dimensions. J. Chem. Phys. 2002, 117, 5780−5784. (21) Gózd́ ź, W.; Ciach, A. Critical point calculation for binary mixtures of symmetric non-additive hard disks. Condens. Matter Phys. 2016, 19, 13002. (22) Kofke, D. A.; Glandt, E. D. Monte Carlo simulation of multicomponent equilibria in a semigrand canonical ensemble. Mol. Phys. 1988, 64, 1105−1131. (23) Metropolis, N.; Ulam, S. M. The Monte Carlo Method. J. Am. Stat. Assoc. 1949, 44, 335−341. (24) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087−1092. (25) Florek, K.; Łukaszewicz, J.; Perkal, J.; Steinhaus, H.; Zubrzycki, S. Sur la liaison et la division des points d’un ensemble fini. Colloquium Math. 1951, 2, 282−285. (26) Florek, K.; Łukaszewicz, J.; Perkal, J.; Steinhaus, H.; Zubrzycki, S. Taksonomia Wrocławska. Przegl. antrop. 1951, 17, 193−211. (27) Sibson, R. SLINK: an optimally efficient algorithm for the single-link cluster method. Computer Journal (British Computer Society) 1973, 16, 30−34. (28) Frenkel, D. Speed-up of Monte Carlo simulations by sampling of rejected states. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 17571− 17575. (29) Hoover, W. G.; Ree, F. H. Melting Transition and Communal Entropy for Hard Spheres. J. Chem. Phys. 1968, 49, 3609−3617. (30) Fisher, M. E.; Barber, M. Scaling Theory for Finite-Size Effects in the Critical Region. Phys. Rev. Lett. 1972, 28, 1516−1519. (31) Domb, C., Lebowitz, J., Eds. Phase Transitions and Critical Phenomena, 8th ed.; Academic Press: London, 1983. (32) Wilding, N. B. Critical-point and coexistence-curve properties of the Lennard-Jones fluid: A finite-size scaling study. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 52, 602−611. (33) Binder, K. Critical Properties from Monte Carlo Coarse Graining and Renormalization. Phys. Rev. Lett. 1981, 47, 693−696. (34) Kamieniarz, G.; Blote, H. W. J. Universal ratio of magnetization moments in two-dimensional Ising models. J. Phys. A: Math. Gen. 1993, 26, 201.

F

DOI: 10.1021/acs.langmuir.7b02228 Langmuir XXXX, XXX, XXX−XXX