Fluid-solid coexistence from two-phase simulations: Binary colloidal mixtures and square well systems G. Arlette Méndez-Maldonado, Gustavo A. Chapela, José Adrián Martínez-González, José Antonio Moreno, Enrique Díaz-Herrera, and José Alejandre Citation: The Journal of Chemical Physics 142, 054501 (2015); doi: 10.1063/1.4906424 View online: http://dx.doi.org/10.1063/1.4906424 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/5?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Note: Effect of fluid phase compositions on the formation of substitutionally ordered solid phases from a binary mixture of oppositely charged colloidal suspensions J. Chem. Phys. 139, 086101 (2013); 10.1063/1.4819952 Determination of favorable inter-particle interactions for formation of substitutionally ordered solid phases from a binary mixture of oppositely charged colloidal suspensions J. Chem. Phys. 138, 174504 (2013); 10.1063/1.4802784 Monte Carlo cluster algorithm for fluid phase transitions in highly size-asymmetrical binary mixtures J. Chem. Phys. 133, 194102 (2010); 10.1063/1.3495996 Statics and dynamics of colloid-polymer mixtures near their critical point of phase separation: A computer simulation study of a continuous Asakura–Oosawa model J. Chem. Phys. 130, 064906 (2009); 10.1063/1.3071197 Grand canonical Monte Carlo simulation of a model colloid–polymer mixture: Coexistence line, critical behavior, and interfacial tension J. Chem. Phys. 121, 3253 (2004); 10.1063/1.1773771
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
THE JOURNAL OF CHEMICAL PHYSICS 142, 054501 (2015)
Fluid-solid coexistence from two-phase simulations: Binary colloidal mixtures and square well systems G. Arlette Méndez-Maldonado,1 Gustavo A. Chapela,2,a) José Adrián Martínez-González,2,b) José Antonio Moreno,2 Enrique Díaz-Herrera,2 and José Alejandre1
1
Departamento de Química, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México Distrito Federal, Mexico 2 Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México Distrito Federal, Mexico
(Received 3 October 2014; accepted 9 January 2015; published online 2 February 2015) Molecular dynamics simulations are performed to clarify the reasons for the disagreement found in a previous publication [G. A. Chapela, F. del Río, and J. Alejandre, J. Chem. Phys. 138(5), 054507 (2013)] regarding the metastability of liquid-vapor coexistence on equimolar charged binary mixtures of fluids interacting with a soft Yukawa potential with κσ = 6. The fluid-solid separation obtained with the two-phase simulation method is found to be in agreement with previous works based on free energy calculations [A. Fortini, A.-P. Hynninen, and M. Dijkstra, J. Chem. Phys. 125, 094502 (2006)] only when the CsCl structure of the solid is used. It is shown that when pressure is increased at constant temperature, the solids are amorphous having different structures, densities, and the diagonal components of the pressure tensor are not equal. A stable low density fluid-solid phase separation is not observed for temperatures above the liquid-vapor critical point. In addition, Monte Carlo and discontinuous molecular dynamics simulations are performed on the square well model of range 1.15σ. A stable fluid-solid transition is observed above the vapor-liquid critical temperature only when the solid has a face centered cubic crystalline structure. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4906424]
I. INTRODUCTION
Since the work by Gast et al.1 on protein induced phase separation in a nonaqueous colloidal suspension, where they predicted that for a short ranged potential of less than a third the size of the hard sphere, no liquid-vapor equilibrium envelope existed. This initial work produced several important works that contributed to the establishment of the nonexistence of the liquid-vapor equilibrium and its critical point or at least to its classification of metastable with respect to the solid-vapor equilibrium envelope. In their computer simulations of dipolar fluids, van Leeuwen and Smit,2 based on this seminal work by Gast et al.,1 affirmed that if the range of the attractive potential is too small, the (metastable) critical temperature is lower than the triple temperature, and hence no vapor-liquid equilibrium is observed. Hagen and Frenkel,3 based on the same work by Gast et al.,1 stated that for short ranged attraction, the liquid phase disappears and only the solid and the fluid phases remained. From their own simulation results, they estimated that the liquid-vapor coexistence curve disappeared when the range of the attractive part of the Yukawa potential is less than approximately one-sixth of the hard-core diameter. These results are further supported by theoretical,4 experimental,5,6 a)Author to whom correspondence should be addressed. Electronic mail:
[email protected]
b)Present address: Institute for Molecular Engineering, The University
of Chicago, Chicago, Illinois 60637, USA. 0021-9606/2015/142(5)/054501/6/$30.00
and simulation work.7 Foffi et al.8 using a hybrid theoretical method determined the equilibrium phase diagram for an attractive Yukawa fluid. They expressed some concern about the minimum limit range of the attractive potential where the liquid-vapor equilibrium becomes metastable. Smaller values of the range of the attractive potential were reported by Mederos and Navascués9 and a larger one by Shukla10 The range found by Shukla doubles the value reported by Hagen and Frenkel3 and Foffi et al.8 For systems with long ranged interactions, a high density fluid-solid transition is found for temperatures above the critical point of the vapor-liquid equilibrium. That is the case of full Lennard-Jones11 fluids, the square well (SW) potential with range of 1.5 molecular diameters,12 the screened Yukawa,13,14 and the restricted primitive model (RPM) for ionic fluids.15,16 Apart from the three phases (vapor, liquid, and solid), the phase diagram contains a well defined vaporliquid critical and triple points. When the interaction range is very short, a stable low density fluid-solid transition, above the vapor-liquid critical point, has been reported. Liu et al.,12 Orkoulas,17 and Pagan and Gunton18 obtained results for the square well with attraction range of 1.15σ, where σ is the hard sphere diameter. Dijkstra13 analyzed the stability of the vapor-liquid transition of a single component of attractive hard core Yukawa fluids while Fortini et al.14 and Hynninen et al.16 that of the binary charged Yukawa colloids. They found that the vapor-liquid equilibrium was metastable with respect to the fluid-solid equilibrium for κσ > 7 and κσ > 4.5 for the one component and binary mixtures, respectively.
142, 054501-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-2
Méndez-Maldonado et al.
For mixtures of attractive and repulsive Yukawas some further debate has taken place between Caballero et al.19 and Fortini et al.14 with respect to the minimum limit range of the potential where the liquid-vapor equilibrium becomes metastable. In a previous work, Chapela et al.20 reviewed several of the reported results for fluids interacting with the Yukawa model. They performed NVT (constant number of particles, volume, and temperature) simulations using continuous molecular dynamics and discontinuous molecular dynamics (DMD) to determine the vapor-liquid phase equilibrium and surface tension for the one component attractive Yukawa fluid with values of κσ from 1.8 to 15 and also for binary mixtures of colloidal particles with interaction parameters κσ from 1.8 to 20. The results of surface tension in both cases were positive for all ranges of interaction. There was not found a signal of metastability. In the same work, Chapela et al.20 performed NPT simulations (constant number of particle, pressure, and temperature) for binary soft Yukawa mixtures for a short ranged interaction model using κσ = 6 in order to analyze the vapor-solid, vapor-liquid, and fluid-solid transitions. From the compression of one phase systems using NPT simulations at a given temperature, the vapor-liquid transition was found, as expected, at densities slightly higher compared to the equilibrium conditions obtained from slab simulations using the NVT ensemble. The density, for temperatures above the vapor-liquid critical point, increased monotonically as pressure increased until a high density fluid-solid transition was found as occurs in systems with long ranged interactions. In those calculations, there was not an evidence of a low density fluid-solid transition, in contrast with results reported by Fortini et al.14 To the best of our knowledge, the stability of the high density fluid-solid transition obtained from NPT simulations has not been discussed for short ranged potentials. According to NPT results at T ∗ = 0.18, obtained by Chapela et al.20 for the binary mixture of colloids with κσ = 6, the fluid density where the fluid-solid transition occurs is 0.8 while the value reported by Fortini et al. is around 0.1. The differences, as it is shown in this work, are because the solid obtained using NPT simulations do not have any crystalline structure. The main goal of this work is to clarify the differences found for the fluid-solid transition from NPT and NVT simulations with results obtained using free energy calculations or other methods to study phase equilibria. The results for the binary colloidal mixtures are compared with the free energy calculations reported by Fortini et al.14 The components of the pressure tensor are used in isotropic volume fluctuations NPT simulations to determine if the solid is mechanically stable where the diagonal components have to be the same.21,22 The agreement between the two-phase method in the NVT ensemble and that of the free energy route is obtained when the solid has a CsCl structure. The same exercise is performed for the SW potential with range 1.15σ using Monte Carlo (MC) in the NPT ensemble and DMD at constant temperature. The NVT simulation results show that the low density fluid-solid equilibrium is stable and above the vapor-liquid critical point only when the solid has a faced center cubic (FCC) crystalline structure. This is the first time
J. Chem. Phys. 142, 054501 (2015)
the phase equilibria for the fluid-solid transition is determined for short ranged interactions using the NVT ensemble. The work is organized as follows. The soft Yukawa potential for binary mixtures and square well model are described in Sec. II together with simulation details. The one phase and two-phases simulation results are discussed in Sec. III. Finally, conclusions and references are given.
II. POTENTIAL MODELS AND SIMULATION DETAILS
The interaction between particles of binary mixture of colloids is given by the soft Yukawa model (SYM) and with the square well potential for those from the one component systems. A. The soft Yukawa potential
The screened soft primitive model23 is used to simulate equimolar mixtures of particles carrying opposite charge. The interaction between two particles α and β at a distance r α β is ) 225 ( qα qβ e−κ(r α β /σ α β −1) σ + Uα β (r α β ) = f min, (1) r α β /σα β r α β where κ −1 is a measure of the range of interaction in dimensionless units and qα is the charge of particle α. All the particles have the same size σ, and all the simulations are performed using κσ = 6. The factor f min = 1.075 is included to have a potential which is around zero at r α β = σ and −1 at the minimum of the attractive interactions. The reduced charges of positive and negative particles have a value of 1.0 and −1.0, respectively. The short-ranged soft model has been used before to study interfacial properties of the restricted primitive model24,25 and ions with asymmetry in size and charge.26 The potential defined in Eq. (1) includes explicitly the f min factor required to have a potential close to zero at σ with a reduced energy of −1.0 for the attractive interactions. In Refs. 23–26 is stated that the parameters were chosen to have both conditions, when the simulations were carried out the f min was included. The equations of motion are solved using the velocity Verlet algorithm with a reduced time step of 0.0005. The Nosé-Hoover chain of thermostats and barostats27,28 is applied to impose the external temperature and pressure. Full periodic boundary conditions are used in all directions. The NPT ensemble is used to carry out one phase simulations to determine the conditions where vapor-fluid and fluid-solid transitions occur. The liquid-vapor and fluid-solid phases in equilibrium are studied using NVT simulations combined with the slab geometry where the high density phase (solid) is surrounded by other phase of a lower value (vapor or liquid). The results are obtained using a parallelized program written in fortran 90. Two initial configurations, FCC and cubic (CsCl), are used in NPT simulations to analyze their effect on the equation of state of the liquid and solid phases. Simulations with 512 and 1024 are performed using the FCC and cubic arrays with CsCl structure, respectively, for the one phase simulations. The potential is truncated at 3.8σ and the density is obtained
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-3
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
for at least 2 × 107 time steps after an equilibration period of half that time. The two-phase simulations, vapor-liquid or fluid-solid equilibria, are performed using the NVT ensemble with a FCC and CsCl crystalline arrays to determine the initial configuration. The high density phase, liquid or solid, is allocated in the middle of the simulation cell surrounded by a vapor or a liquid. The fluid-solid simulations are performed starting with a solid having a CsCl structure. The systems contain 2048 particles in non-cubic simulation cells with typical dimensions of L x = 8σ and L z = 3L x . The truncation distance is also 3.8σ. The average properties are obtained for at least 3 × 108 time steps after an equilibration period of 2 × 107 time steps. B. The square well potential
The interaction potential for the SW model with a range of interaction λ is −ϵ, σ u(r) = 0,
λ, λ,
and Hynninen et al.16 is used as a guide to choose the total density and temperature for the vapor-liquid and fluid-solid using the two-phase simulations method. The MC results in the NPT ensemble using the SW potential are obtained for the vapor-liquid and fluid-solid at constant pressure and several temperatures. The DMD simulation results using the NVT ensemble are also given for the fluid-solid transition. A. Binary colloidal mixtures
In a previous work, Chapela et al.20 reported the vaporliquid phase diagram for the SYM with interaction parameter κσ = 6, and the critical temperature was estimated to be 0.175. NVT simulations were performed in that work with a liquid slab surrounded by vapor. The average densities obtained with that method correspond to equilibrium conditions. Chapela et al.20 also reported a fluid-solid like transition using NPT simulations.
(2)
where r is the distance between particles, and σ and ϵ are the diameter of the hard core and the depth of the square well. One phase simulations are performed using the MC method in the NPT ensemble29 with cubic boxes where the fluctuations in the volume of the simulation box are implemented by changing, simultaneously, all the sides of the box by the same amount. The NPT simulations29,30 are performed with N particles to locate the vapor-liquid and fluid-solid transitions at different temperatures. It is known that the jump in pressure does not occur at the orthobaric curve but instead at some intermediate point between it and the spinodal curve. DMD in the canonical ensemble20,31–34 is used to perform the two-phase simulations, vapor-liquid or fluid-solid. A parallelepiped box is used with L x = L y and L z /L x = 2 or 3. Periodic boundary conditions are used in all directions. Simulations are run in blocks of 25 × 106 collisions. The minimum number of blocks performed is 79 and the maximum is 575. Most runs are longer than 150 blocks. The analysis of the properties given in this work is done with, at least, half of the collisions performed. Reduced units are used in this work: distance r ∗ = r/σ, energy u∗ = u/ϵ (where ϵ = U(r min)), temperature T ∗ = kT/ϵ, density ρ∗ = ρσ 3, time ∆t ∗ = ∆t(ϵ/mσ 2)1/2, pressure p∗ = pσ 3/ϵ, λ∗ = λ/σ. The reduced masses of the particles used for the DMD simulations of the square well model are 1.
1. SYM one phase simulations
NPT simulations are carried out in this work to analyze the solid phase with more detail and to obtain the density as a function of pressure for temperatures ranging from T ∗ = 0.12 to 0.185. The range of temperature covers the vapor-liquid and fluid-solid phase transitions. The initial configuration started with particles randomly distributed in a FCC crystalline array. The results are shown in Figure 1. At constant temperature, as the pressure increases, the density also increases. The vapor-solid transition is observed for T ∗ = 0.12 and 0.14 while the vapor-liquid phase separation is found for temperatures between 0.15 and 0.17. The triple point temperature lies between T ∗ = 0.14 and 0.15. The high density fluid-solid phase transition is found for temperatures starting at 0.15. A solidfluid transition is observed below the triple point. It can be seen from Fig. 1 that at temperatures above the vapor-liquid critical point, the density increases monotonically with pressure and there is a jump in the high density region, resembling a first order phase transition. Amorphous solids with different
III. RESULTS
This section contains results for one phase and two-phase simulations using the NPT and NVT ensembles. The stability of the solid phase for the binary colloidal mixtures is analyzed by the use of the diagonal components of the pressure tensor of solids having amorphous and CsCl structures. The amorphous solid is obtained by compressing the liquid at constant temperature until a fluid-solid like transition is observed. A direct comparison between two-phase simulations and free energy calculations allows showing that the results for the fluid-solid equilibrium from both methods are in excellent agreement. The phase diagram reported by Chapela et al.20
FIG. 1. Density as a function of pressure obtained in one phase NPT simulations at temperatures shown in the figure for κσ = 6. Open and filled symbols are for initial configurations from a FCC and CsCl crystalline arrays, respectively.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-4
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
TABLE I. Diagonal components of the pressure tensor for solids obtained from NPT simulations at different temperatures at the end of the compression procedure starting from a vapor phase. The initial configuration is chosen to be a solid with a FCC or CsCl crystalline structure. Initial array
T∗
FCC FCC FCC FCC FCC FCC CsCl CsCl
0.120 0.150 0.170 0.175 0.180 0.185 0.160 0.185
ρ∗
Px x
Py y
Pz z
0.948 0.873 1.005 1.072 1.084 1.054 1.040 1.062
0.33 0.13 −0.31 0.29 1.10 0.45 0.39 0.49
−0.14 0.17 2.02 1.48 0.38 0.17 0.24 0.73
−0.19 0.004 0.69 3.32 0.33 1.48 0.87 0.27
densities are observed. Additional NPT simulations at T ∗ = 0.16 and T ∗ = 0.185 are performed using a CsCl crystalline array as initial configuration to determine the effect of the initial configuration on the stability of the solids. The results are also shown in Fig. 1 and the solids are the same as those obtained using a FCC crystalline array as initial configuration. To analyze the stability of the solids, the diagonal components of the pressure tensor are calculated. In a stable solid, the three components have to be the same.21 The results given in Table I show that, at the same temperature, the three components are very different, and therefore, the solids are not mechanically stable. The diagonal components of the pressure tensor for systems starting with a CsCl initial configuration are also different, see Table I. In NPT simulations, the vapor-liquid and fluid-solid transitions occur at conditions different from equilibrium; hence, in the compression, the vapor and liquid densities are slightly greater than values at equilibrium. The Fig. 2 shows the NPT results for the vapor-liquid and fluid-solid transitions. It is noted that the values of the density plotted in Fig. 2 are those obtained where the transitions occur. The vapor-liquid phase
FIG. 2. Density as a function of temperature for Yukawa mixtures with κσ = 6. The liquid-vapor coexistence obtained by Chapela et al.20 (open squares) and from this work (filled diamonds) using NVT and NPT ensembles, respectively, is shown. The critical point is shown by a star. The open circles are results for the unstable fluid-solid coexistence Yukawa binary mixtures with κσ = 6 from NPT simulations. The RPM results for the fluids-solid equilibrium are taken from Ref. 15.
diagram obtained by Chapela et al.20 using NVT simulations at equilibrium is also shown. The results from this work and those from Chapela et al.,20 as expected, are slightly different. We did not attempt to obtain the equilibrium densities and pressures from NPT simulations in the high density fluid-solid transition due to the difference in the components of the pressure tensor of the solid phase. The densities at different temperatures for the unstable fluid-solid transition for the binary mixtures of colloid particles which interact with the soft Yukawa potential are shown in Fig. 2, and they are compared with those of the stable equilibrium found in ionic fluids described by the RPM.15 Although the density values from NPT simulations are closed, the soft Yukawa mixture with κσ = 6 and RPM systems have a completely different fluid-solid equilibrium. 2. SYM two-phase simulations
The equilibrium of the fluid-solid transition above the vapor-liquid critical temperature is also analyzed using NVT simulations performed in a slab geometry, in which, the CsCl crystal structure is surrounded by a fluid of lower density. Several simulations are performed at T ∗ = 0.185 and total densities from 0.43 to 0.8 using 2048 ions. The density profiles show that the fluid-solid equilibrium is stable after 380 × 106 time steps for all the cases analyzed in this work, see Fig. 3. All the simulations give the same vapor and solid densities, a greater amount of vapor is observed when the length of the simulation box in the Z direction is increased in order to decrease the total density of the system. The fluid-solid (CsCl) phase separation is observed and it is in equilibrium for conditions above the vapor-liquid critical point. The twophase simulation method, where a fluid coexists with a CsCl solid, is used to obtain the coexisting densities for higher temperatures, a typical density profile at T ∗ = 0.35 is shown in Fig. 4 for a total density of 0.95, where the fluid and solid regions are well defined after the system has evolved for 50 × 106 steps. The average of coexisting densities is obtained from the density profiles on the bulk regions. The components of the pressure tensor are calculated at the end of the NVT simulations, the results are given in
FIG. 3. Density profiles for a Yukawa colloidal mixture with κσ = 6 at T∗ = 0.185 using as initial configuration a disorder fluid at low density and a solid with a CsCl structure. The colors are used to indicate the total density: magenta (ρ ∗ = 0.8), blue (ρ ∗ = 0.68), green (ρ ∗ = 0.51), red (ρ ∗ = 0.48), and black (ρ ∗ = 0.43).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-5
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
FIG. 4. Density profile for a Yukawa colloidal binary mixture with κσ = 6 at T∗ = 0.3 using a solid with a CsCl structure as initial configuration. A slab of solid is surrounded by two regions of fluid with a total density of 0.95.
Table II. The Px x and Py y are equal, as it should be in a well equilibrated interface. Finally, the fluid-solid phase diagram of the binary Yukawa mixtures with κσ = 6 is obtained using NVT simulations where a slab of CsCl solid is surrounded by a fluid with different densities depending on the temperature. The results from this work for the fluid-solid equilibrium are in excellent agreement with those reported by Fortini et al.14 in the fluid branch for temperatures up to the vapor-liquid critical point and by Hynninen et al.16 in both branches for temperatures above 0.2. The solid branch is at a lower density than the reported by both works.14,16 See Fig. 5 for a schematic representation of these results. In addition, the vapor-liquid phase diagram obtained by Chapela et al.20 and Fortini et al.14 is also shown. B. Square well potential
Extensive simulations are performed to analyze the vaporliquid and fluid-solid transitions for the SW model with λ = 1.15σ. One phase MC simulations in the NPT ensemble with isotropic fluctuation of volume are performed for nine isobars. Two-phase simulations using the DMD on the NVT ensemble are also carried out to obtain the low density fluidsolid equilibrium using the FCC structure for the solid phase.
FIG. 5. Phase diagram of Yukawa mixtures with κσ = 6. The vapor-liquid coexistence from Chapela et al.20 and Fortini et al.14 is shown with filled squares and open circles, respectively. The fluid-solid coexistence results are also shown: Hynninen et al.16 (filled triangles), Fortini et al.14 (filled circles), and this work using NVT simulations (filled squares).
pressures 0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.5. If the initial conditions of the simulation are within the spinodal region, a violent separation of the liquid and the vapor phases takes place. The vapor-liquid and high density fluid-solid transitions are observed at the lowest pressure. For reduced pressures above or equal 0.5, only the high density-solid transition is observed. Solids with different densities and structures, as in the case of binary Yukawa mixtures, are found. Up to this point, the results look similar to those of long ranged interactions where a low density fluid-solid transition is not observed. The existence of that transition published by several authors12,17,18 is not obtained using NPT simulations. The stability of the solid phase is determined by the pressure tensor components, for example, at T ∗ = 0.58, their values are Py y = 0.226 and Pz z = 0.099. As in the case of the binary colloidal mixtures, they were not the same.
1. SW one phase simulations
The MC results at constant pressure and different temperatures are shown in Fig. 6 for isobars with reduced TABLE II. Components of the pressure tensor for the fluid-solid interfaces obtained from NVT simulations at different temperatures. The total density of the system is also given. T∗ 0.185 0.185 0.185 0.200 0.300 0.400
ρ∗
Px x
Py y
Pz z
0.43 0.68 0.80 0.80 0.80 0.95
0.075 0.122 0.145 0.250 1.020 2.480
0.075 0.122 0.145 0.250 1.020 2.480
0.006 0.006 0.006 0.050 0.800 1.680
FIG. 6. MC NPT results for the SW potential. The fluid-solid and vaporliquid equilibrium envelopes are shown. From bottom to top, the pressures are 0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.5. The filled circles are vapor-liquid curves from this work and open circles correspond to Orkoulas.17 The dashed curve represents the place where the liquid jumps to the solid phase. The low density fluid-solid equilibrium is not observed from NPT simulations.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-6
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
solids have a CsCl structure, the pressure components are the same. The procedure to obtain the fluid-solid phase diagram should include the analysis of the components of the pressure tensor of the solid phase to establish a mechanical equilibrium necessary condition in addition to the chemical potential, pressure, and temperature. The use of an ideal crystal in the thermodynamic integration method to draw the fluid-solid boundaries leads to a solid with equal diagonal components of the pressure tensor. ACKNOWLEDGMENTS
FIG. 7. Results for the SW using DMD method in the NVT ensemble. Filled squares represent vapor-solid equilibrium curves started with a FCC crystal obtained from runs at densities of 0.641 and 0.9 obtained in this work. Open circles are taken from Orkoulas17 and open triangles from Liu et al.12 The vapor-liquid points obtained in this work are shown with filled circles.
2. SW two-phase simulations
To test the stability of the vapor-liquid phase equilibrium for short ranged interactions, a perturbation from the solution at high density is carried out. If the state is stable, any perturbation must reach the same equilibrium point. The maximum perturbation consists in starting the NVT DMD simulation runs from a solid, obtained from MC NPT simulations. The solid is surrounded by vacuum to obtain a two-phase system. As the system evolves in time, the solid melts and a homogeneous fluid is obtained at the end of the simulation. The reason why the fluid-solid is not formed is because the solid is not stable, as explained above. The crystal structure for the SW model is FCC. When a FCC crystal is surrounded by a low or high density fluid, a low density-crystal equilibrium is obtained. The phase equilibrium results are shown in Fig. 7. The results agree well, within simulation error, with those published by Liu et al.12 and Orkoulas17 using other simulation methods. Small differences at high temperatures are found in the fluid phase where the densities from this work are slightly larger. IV. CONCLUSIONS
Results from this work for short ranged interactions show that the low density fluid-solid equilibrium is only found when the solid has a crystalline structure, CsCl for the binary colloidal mixtures with κσ = 6 and FCC for the square well potential with λ = 1.15σ. The results for binary colloidal mixture using the two-phase simulation method in the NVT ensemble give the same values for the fluid-solid equilibrium as those reported using the free energy calculations14,16 only when the solid phase has a CsCl structure. For the SW fluid, DMD results using the NVT ensemble were around the same as those published by Liu et al.12 and Orkoulas17 using other simulation methods only when the solid has a FCC structure. The solids obtained in NPT simulations at all temperatures have different diagonal components of the pressure tensor, and therefore, they are not mechanically at equilibrium. Therefore, the high density fluid-solid transition is not stable. When the
The authors thank to Laboratorio de Supercómputo for allocation time. J.A. thanks Enrique Velasco for providing us the program to build the CsCl solid structure. We thank the referees for helpful suggestions. 1A.
P. Gast, C. K. Hall, and W. B. Russel, J. Colloid Interface Sci. 96, 251 (1983). 2M. E. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 3991 (1993). 3M. H. J. Hagen and D. Frenkel, J. Chem. Phys. 101(5), 4093 (1994). 4H. N. W. Lekkerkerker, W. C. K. Poon, P. N. Pusey, A. Stroobants, and P. B. Warn, Europhys. Lett. 20, 559 (1992). 5P. N. Pusey, W. C. K. Poon, S. M. Ilett, and P. Barlett, J. Phys.: Condens. Matter 6, A29 (1994). 6E. Leal Calderon, J. Bibette, and J. Biais, Europhys. Lett. 23, 653 (1993). 7E. J. Meijer and D. Frenkel, Phys. Rev. Lett. 67, 1110 (1991). 8G. Foffi, G, D. McCullagh, A. Lawlor, E. Zaccarelli, K. A. Dawson, F. Sciortino, P. Tartaglia, D. Pini, and G. Stell, Phys. Rev. E 65, 031407 (2002). 9L. Mederos and G. Navascués, J. Chem. Phys. 101, 9841 (1994). 10K. P. Shukla, J. Chem. Phys. 112, 10358 (2000). 11M. A. Barroso and A. L. Ferreira, J. Chem. Phys. 116, 7145 (2002). 12H. Liu, S. Garde, and S. Kumara, J. Chem. Phys. 123, 174505 (2005). 13M. Dijkstra, Phys. Rev. E 66, 021402 (2002). 14A. Fortini, A.-P. Hynninen, and M. Dijkstra, J. Chem. Phys. 125, 094502 (2006). 15C. Vega, J. L. F. Abascal, C. McBride, and F. Bresme, J. Chem. Phys. 119, 964 (2003). 16A.-P. Hynninen, M. E. Leunnissen, A. van Blaaderen, and M. Dijkstra, Phys. Rev. Lett. 96, 018303 (2006). 17G. Orkoulas, J. Chem. Phys. 133, 111104 (2010). 18D. L. Pagan and J. D. Gunton, J. Chem. Phys. 122, 184515 (2005). 19J. B. Caballero, A. M. Puertas, A. Fernández-Barbero, F. J. de las Nieves, J. M. Romero-Enrique, and L. F. Rull, J. Chem. Phys. 124, 054909 (2006). 20G. A. Chapela, F. del Río, and J. Alejandre, J. Chem. Phys. 138(5), 054507 (2013). 21M. P. Allen and A. J. Master, Mol. Phys. 79, 277 (1993). 22H. Domínguez, E. Velasco, and J. Alejandre, Mol. Phys. 100, 2739 (2002). 23F. Bresme, M. González-Melchor, and J. Alejandre, J. Phys.: Condens. Matter. 17, S3301 (2005). 24M. González-Melchor, J. Alejandre, and F. Bresme, Phys. Rev. Lett. 90, 135506 (2003). 25M. González-Melchor, F. Bresme, and J. Alejandre, J. Chem. Phys. 122, 104710 (2005). 26J. Alejandre, F. Bresme, and M. González-Melchor, Mol. Phys. 107, 357 (2009). 27M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, J. Phys. A: Math. Gen. 39, 5629 (2006). 28M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulations (Oxford University Press, New York, 2010). 29D. Frenkel and B. Smith, Understanding Molecular Simulation from Algorithms to Applications, 2nd ed. (Academic Press, San Diego, 1996). 30P. Gregory and J. C. Schug, Mol. Phys. 82(4), 677 (1994). 31B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, 459 (1959). 32D. C. Rapaport, J. Chem. Phys. 71, 3299 (1979). 33G. A. Chapela, S. E. Martínez-Casas, and J. Alejandre, Mol. Phys. 53, 139 (1984). 34P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford University Press, New York, 1987).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
THE JOURNAL OF CHEMICAL PHYSICS 142, 054501 (2015)
Fluid-solid coexistence from two-phase simulations: Binary colloidal mixtures and square well systems G. Arlette Méndez-Maldonado,1 Gustavo A. Chapela,2,a) José Adrián Martínez-González,2,b) José Antonio Moreno,2 Enrique Díaz-Herrera,2 and José Alejandre1
1
Departamento de Química, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México Distrito Federal, Mexico 2 Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, 09340 México Distrito Federal, Mexico
(Received 3 October 2014; accepted 9 January 2015; published online 2 February 2015) Molecular dynamics simulations are performed to clarify the reasons for the disagreement found in a previous publication [G. A. Chapela, F. del Río, and J. Alejandre, J. Chem. Phys. 138(5), 054507 (2013)] regarding the metastability of liquid-vapor coexistence on equimolar charged binary mixtures of fluids interacting with a soft Yukawa potential with κσ = 6. The fluid-solid separation obtained with the two-phase simulation method is found to be in agreement with previous works based on free energy calculations [A. Fortini, A.-P. Hynninen, and M. Dijkstra, J. Chem. Phys. 125, 094502 (2006)] only when the CsCl structure of the solid is used. It is shown that when pressure is increased at constant temperature, the solids are amorphous having different structures, densities, and the diagonal components of the pressure tensor are not equal. A stable low density fluid-solid phase separation is not observed for temperatures above the liquid-vapor critical point. In addition, Monte Carlo and discontinuous molecular dynamics simulations are performed on the square well model of range 1.15σ. A stable fluid-solid transition is observed above the vapor-liquid critical temperature only when the solid has a face centered cubic crystalline structure. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4906424]
I. INTRODUCTION
Since the work by Gast et al.1 on protein induced phase separation in a nonaqueous colloidal suspension, where they predicted that for a short ranged potential of less than a third the size of the hard sphere, no liquid-vapor equilibrium envelope existed. This initial work produced several important works that contributed to the establishment of the nonexistence of the liquid-vapor equilibrium and its critical point or at least to its classification of metastable with respect to the solid-vapor equilibrium envelope. In their computer simulations of dipolar fluids, van Leeuwen and Smit,2 based on this seminal work by Gast et al.,1 affirmed that if the range of the attractive potential is too small, the (metastable) critical temperature is lower than the triple temperature, and hence no vapor-liquid equilibrium is observed. Hagen and Frenkel,3 based on the same work by Gast et al.,1 stated that for short ranged attraction, the liquid phase disappears and only the solid and the fluid phases remained. From their own simulation results, they estimated that the liquid-vapor coexistence curve disappeared when the range of the attractive part of the Yukawa potential is less than approximately one-sixth of the hard-core diameter. These results are further supported by theoretical,4 experimental,5,6 a)Author to whom correspondence should be addressed. Electronic mail:
[email protected]
b)Present address: Institute for Molecular Engineering, The University
of Chicago, Chicago, Illinois 60637, USA. 0021-9606/2015/142(5)/054501/6/$30.00
and simulation work.7 Foffi et al.8 using a hybrid theoretical method determined the equilibrium phase diagram for an attractive Yukawa fluid. They expressed some concern about the minimum limit range of the attractive potential where the liquid-vapor equilibrium becomes metastable. Smaller values of the range of the attractive potential were reported by Mederos and Navascués9 and a larger one by Shukla10 The range found by Shukla doubles the value reported by Hagen and Frenkel3 and Foffi et al.8 For systems with long ranged interactions, a high density fluid-solid transition is found for temperatures above the critical point of the vapor-liquid equilibrium. That is the case of full Lennard-Jones11 fluids, the square well (SW) potential with range of 1.5 molecular diameters,12 the screened Yukawa,13,14 and the restricted primitive model (RPM) for ionic fluids.15,16 Apart from the three phases (vapor, liquid, and solid), the phase diagram contains a well defined vaporliquid critical and triple points. When the interaction range is very short, a stable low density fluid-solid transition, above the vapor-liquid critical point, has been reported. Liu et al.,12 Orkoulas,17 and Pagan and Gunton18 obtained results for the square well with attraction range of 1.15σ, where σ is the hard sphere diameter. Dijkstra13 analyzed the stability of the vapor-liquid transition of a single component of attractive hard core Yukawa fluids while Fortini et al.14 and Hynninen et al.16 that of the binary charged Yukawa colloids. They found that the vapor-liquid equilibrium was metastable with respect to the fluid-solid equilibrium for κσ > 7 and κσ > 4.5 for the one component and binary mixtures, respectively.
142, 054501-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-2
Méndez-Maldonado et al.
For mixtures of attractive and repulsive Yukawas some further debate has taken place between Caballero et al.19 and Fortini et al.14 with respect to the minimum limit range of the potential where the liquid-vapor equilibrium becomes metastable. In a previous work, Chapela et al.20 reviewed several of the reported results for fluids interacting with the Yukawa model. They performed NVT (constant number of particles, volume, and temperature) simulations using continuous molecular dynamics and discontinuous molecular dynamics (DMD) to determine the vapor-liquid phase equilibrium and surface tension for the one component attractive Yukawa fluid with values of κσ from 1.8 to 15 and also for binary mixtures of colloidal particles with interaction parameters κσ from 1.8 to 20. The results of surface tension in both cases were positive for all ranges of interaction. There was not found a signal of metastability. In the same work, Chapela et al.20 performed NPT simulations (constant number of particle, pressure, and temperature) for binary soft Yukawa mixtures for a short ranged interaction model using κσ = 6 in order to analyze the vapor-solid, vapor-liquid, and fluid-solid transitions. From the compression of one phase systems using NPT simulations at a given temperature, the vapor-liquid transition was found, as expected, at densities slightly higher compared to the equilibrium conditions obtained from slab simulations using the NVT ensemble. The density, for temperatures above the vapor-liquid critical point, increased monotonically as pressure increased until a high density fluid-solid transition was found as occurs in systems with long ranged interactions. In those calculations, there was not an evidence of a low density fluid-solid transition, in contrast with results reported by Fortini et al.14 To the best of our knowledge, the stability of the high density fluid-solid transition obtained from NPT simulations has not been discussed for short ranged potentials. According to NPT results at T ∗ = 0.18, obtained by Chapela et al.20 for the binary mixture of colloids with κσ = 6, the fluid density where the fluid-solid transition occurs is 0.8 while the value reported by Fortini et al. is around 0.1. The differences, as it is shown in this work, are because the solid obtained using NPT simulations do not have any crystalline structure. The main goal of this work is to clarify the differences found for the fluid-solid transition from NPT and NVT simulations with results obtained using free energy calculations or other methods to study phase equilibria. The results for the binary colloidal mixtures are compared with the free energy calculations reported by Fortini et al.14 The components of the pressure tensor are used in isotropic volume fluctuations NPT simulations to determine if the solid is mechanically stable where the diagonal components have to be the same.21,22 The agreement between the two-phase method in the NVT ensemble and that of the free energy route is obtained when the solid has a CsCl structure. The same exercise is performed for the SW potential with range 1.15σ using Monte Carlo (MC) in the NPT ensemble and DMD at constant temperature. The NVT simulation results show that the low density fluid-solid equilibrium is stable and above the vapor-liquid critical point only when the solid has a faced center cubic (FCC) crystalline structure. This is the first time
J. Chem. Phys. 142, 054501 (2015)
the phase equilibria for the fluid-solid transition is determined for short ranged interactions using the NVT ensemble. The work is organized as follows. The soft Yukawa potential for binary mixtures and square well model are described in Sec. II together with simulation details. The one phase and two-phases simulation results are discussed in Sec. III. Finally, conclusions and references are given.
II. POTENTIAL MODELS AND SIMULATION DETAILS
The interaction between particles of binary mixture of colloids is given by the soft Yukawa model (SYM) and with the square well potential for those from the one component systems. A. The soft Yukawa potential
The screened soft primitive model23 is used to simulate equimolar mixtures of particles carrying opposite charge. The interaction between two particles α and β at a distance r α β is ) 225 ( qα qβ e−κ(r α β /σ α β −1) σ + Uα β (r α β ) = f min, (1) r α β /σα β r α β where κ −1 is a measure of the range of interaction in dimensionless units and qα is the charge of particle α. All the particles have the same size σ, and all the simulations are performed using κσ = 6. The factor f min = 1.075 is included to have a potential which is around zero at r α β = σ and −1 at the minimum of the attractive interactions. The reduced charges of positive and negative particles have a value of 1.0 and −1.0, respectively. The short-ranged soft model has been used before to study interfacial properties of the restricted primitive model24,25 and ions with asymmetry in size and charge.26 The potential defined in Eq. (1) includes explicitly the f min factor required to have a potential close to zero at σ with a reduced energy of −1.0 for the attractive interactions. In Refs. 23–26 is stated that the parameters were chosen to have both conditions, when the simulations were carried out the f min was included. The equations of motion are solved using the velocity Verlet algorithm with a reduced time step of 0.0005. The Nosé-Hoover chain of thermostats and barostats27,28 is applied to impose the external temperature and pressure. Full periodic boundary conditions are used in all directions. The NPT ensemble is used to carry out one phase simulations to determine the conditions where vapor-fluid and fluid-solid transitions occur. The liquid-vapor and fluid-solid phases in equilibrium are studied using NVT simulations combined with the slab geometry where the high density phase (solid) is surrounded by other phase of a lower value (vapor or liquid). The results are obtained using a parallelized program written in fortran 90. Two initial configurations, FCC and cubic (CsCl), are used in NPT simulations to analyze their effect on the equation of state of the liquid and solid phases. Simulations with 512 and 1024 are performed using the FCC and cubic arrays with CsCl structure, respectively, for the one phase simulations. The potential is truncated at 3.8σ and the density is obtained
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-3
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
for at least 2 × 107 time steps after an equilibration period of half that time. The two-phase simulations, vapor-liquid or fluid-solid equilibria, are performed using the NVT ensemble with a FCC and CsCl crystalline arrays to determine the initial configuration. The high density phase, liquid or solid, is allocated in the middle of the simulation cell surrounded by a vapor or a liquid. The fluid-solid simulations are performed starting with a solid having a CsCl structure. The systems contain 2048 particles in non-cubic simulation cells with typical dimensions of L x = 8σ and L z = 3L x . The truncation distance is also 3.8σ. The average properties are obtained for at least 3 × 108 time steps after an equilibration period of 2 × 107 time steps. B. The square well potential
The interaction potential for the SW model with a range of interaction λ is −ϵ, σ u(r) = 0,
λ, λ,
and Hynninen et al.16 is used as a guide to choose the total density and temperature for the vapor-liquid and fluid-solid using the two-phase simulations method. The MC results in the NPT ensemble using the SW potential are obtained for the vapor-liquid and fluid-solid at constant pressure and several temperatures. The DMD simulation results using the NVT ensemble are also given for the fluid-solid transition. A. Binary colloidal mixtures
In a previous work, Chapela et al.20 reported the vaporliquid phase diagram for the SYM with interaction parameter κσ = 6, and the critical temperature was estimated to be 0.175. NVT simulations were performed in that work with a liquid slab surrounded by vapor. The average densities obtained with that method correspond to equilibrium conditions. Chapela et al.20 also reported a fluid-solid like transition using NPT simulations.
(2)
where r is the distance between particles, and σ and ϵ are the diameter of the hard core and the depth of the square well. One phase simulations are performed using the MC method in the NPT ensemble29 with cubic boxes where the fluctuations in the volume of the simulation box are implemented by changing, simultaneously, all the sides of the box by the same amount. The NPT simulations29,30 are performed with N particles to locate the vapor-liquid and fluid-solid transitions at different temperatures. It is known that the jump in pressure does not occur at the orthobaric curve but instead at some intermediate point between it and the spinodal curve. DMD in the canonical ensemble20,31–34 is used to perform the two-phase simulations, vapor-liquid or fluid-solid. A parallelepiped box is used with L x = L y and L z /L x = 2 or 3. Periodic boundary conditions are used in all directions. Simulations are run in blocks of 25 × 106 collisions. The minimum number of blocks performed is 79 and the maximum is 575. Most runs are longer than 150 blocks. The analysis of the properties given in this work is done with, at least, half of the collisions performed. Reduced units are used in this work: distance r ∗ = r/σ, energy u∗ = u/ϵ (where ϵ = U(r min)), temperature T ∗ = kT/ϵ, density ρ∗ = ρσ 3, time ∆t ∗ = ∆t(ϵ/mσ 2)1/2, pressure p∗ = pσ 3/ϵ, λ∗ = λ/σ. The reduced masses of the particles used for the DMD simulations of the square well model are 1.
1. SYM one phase simulations
NPT simulations are carried out in this work to analyze the solid phase with more detail and to obtain the density as a function of pressure for temperatures ranging from T ∗ = 0.12 to 0.185. The range of temperature covers the vapor-liquid and fluid-solid phase transitions. The initial configuration started with particles randomly distributed in a FCC crystalline array. The results are shown in Figure 1. At constant temperature, as the pressure increases, the density also increases. The vapor-solid transition is observed for T ∗ = 0.12 and 0.14 while the vapor-liquid phase separation is found for temperatures between 0.15 and 0.17. The triple point temperature lies between T ∗ = 0.14 and 0.15. The high density fluid-solid phase transition is found for temperatures starting at 0.15. A solidfluid transition is observed below the triple point. It can be seen from Fig. 1 that at temperatures above the vapor-liquid critical point, the density increases monotonically with pressure and there is a jump in the high density region, resembling a first order phase transition. Amorphous solids with different
III. RESULTS
This section contains results for one phase and two-phase simulations using the NPT and NVT ensembles. The stability of the solid phase for the binary colloidal mixtures is analyzed by the use of the diagonal components of the pressure tensor of solids having amorphous and CsCl structures. The amorphous solid is obtained by compressing the liquid at constant temperature until a fluid-solid like transition is observed. A direct comparison between two-phase simulations and free energy calculations allows showing that the results for the fluid-solid equilibrium from both methods are in excellent agreement. The phase diagram reported by Chapela et al.20
FIG. 1. Density as a function of pressure obtained in one phase NPT simulations at temperatures shown in the figure for κσ = 6. Open and filled symbols are for initial configurations from a FCC and CsCl crystalline arrays, respectively.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-4
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
TABLE I. Diagonal components of the pressure tensor for solids obtained from NPT simulations at different temperatures at the end of the compression procedure starting from a vapor phase. The initial configuration is chosen to be a solid with a FCC or CsCl crystalline structure. Initial array
T∗
FCC FCC FCC FCC FCC FCC CsCl CsCl
0.120 0.150 0.170 0.175 0.180 0.185 0.160 0.185
ρ∗
Px x
Py y
Pz z
0.948 0.873 1.005 1.072 1.084 1.054 1.040 1.062
0.33 0.13 −0.31 0.29 1.10 0.45 0.39 0.49
−0.14 0.17 2.02 1.48 0.38 0.17 0.24 0.73
−0.19 0.004 0.69 3.32 0.33 1.48 0.87 0.27
densities are observed. Additional NPT simulations at T ∗ = 0.16 and T ∗ = 0.185 are performed using a CsCl crystalline array as initial configuration to determine the effect of the initial configuration on the stability of the solids. The results are also shown in Fig. 1 and the solids are the same as those obtained using a FCC crystalline array as initial configuration. To analyze the stability of the solids, the diagonal components of the pressure tensor are calculated. In a stable solid, the three components have to be the same.21 The results given in Table I show that, at the same temperature, the three components are very different, and therefore, the solids are not mechanically stable. The diagonal components of the pressure tensor for systems starting with a CsCl initial configuration are also different, see Table I. In NPT simulations, the vapor-liquid and fluid-solid transitions occur at conditions different from equilibrium; hence, in the compression, the vapor and liquid densities are slightly greater than values at equilibrium. The Fig. 2 shows the NPT results for the vapor-liquid and fluid-solid transitions. It is noted that the values of the density plotted in Fig. 2 are those obtained where the transitions occur. The vapor-liquid phase
FIG. 2. Density as a function of temperature for Yukawa mixtures with κσ = 6. The liquid-vapor coexistence obtained by Chapela et al.20 (open squares) and from this work (filled diamonds) using NVT and NPT ensembles, respectively, is shown. The critical point is shown by a star. The open circles are results for the unstable fluid-solid coexistence Yukawa binary mixtures with κσ = 6 from NPT simulations. The RPM results for the fluids-solid equilibrium are taken from Ref. 15.
diagram obtained by Chapela et al.20 using NVT simulations at equilibrium is also shown. The results from this work and those from Chapela et al.,20 as expected, are slightly different. We did not attempt to obtain the equilibrium densities and pressures from NPT simulations in the high density fluid-solid transition due to the difference in the components of the pressure tensor of the solid phase. The densities at different temperatures for the unstable fluid-solid transition for the binary mixtures of colloid particles which interact with the soft Yukawa potential are shown in Fig. 2, and they are compared with those of the stable equilibrium found in ionic fluids described by the RPM.15 Although the density values from NPT simulations are closed, the soft Yukawa mixture with κσ = 6 and RPM systems have a completely different fluid-solid equilibrium. 2. SYM two-phase simulations
The equilibrium of the fluid-solid transition above the vapor-liquid critical temperature is also analyzed using NVT simulations performed in a slab geometry, in which, the CsCl crystal structure is surrounded by a fluid of lower density. Several simulations are performed at T ∗ = 0.185 and total densities from 0.43 to 0.8 using 2048 ions. The density profiles show that the fluid-solid equilibrium is stable after 380 × 106 time steps for all the cases analyzed in this work, see Fig. 3. All the simulations give the same vapor and solid densities, a greater amount of vapor is observed when the length of the simulation box in the Z direction is increased in order to decrease the total density of the system. The fluid-solid (CsCl) phase separation is observed and it is in equilibrium for conditions above the vapor-liquid critical point. The twophase simulation method, where a fluid coexists with a CsCl solid, is used to obtain the coexisting densities for higher temperatures, a typical density profile at T ∗ = 0.35 is shown in Fig. 4 for a total density of 0.95, where the fluid and solid regions are well defined after the system has evolved for 50 × 106 steps. The average of coexisting densities is obtained from the density profiles on the bulk regions. The components of the pressure tensor are calculated at the end of the NVT simulations, the results are given in
FIG. 3. Density profiles for a Yukawa colloidal mixture with κσ = 6 at T∗ = 0.185 using as initial configuration a disorder fluid at low density and a solid with a CsCl structure. The colors are used to indicate the total density: magenta (ρ ∗ = 0.8), blue (ρ ∗ = 0.68), green (ρ ∗ = 0.51), red (ρ ∗ = 0.48), and black (ρ ∗ = 0.43).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-5
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
FIG. 4. Density profile for a Yukawa colloidal binary mixture with κσ = 6 at T∗ = 0.3 using a solid with a CsCl structure as initial configuration. A slab of solid is surrounded by two regions of fluid with a total density of 0.95.
Table II. The Px x and Py y are equal, as it should be in a well equilibrated interface. Finally, the fluid-solid phase diagram of the binary Yukawa mixtures with κσ = 6 is obtained using NVT simulations where a slab of CsCl solid is surrounded by a fluid with different densities depending on the temperature. The results from this work for the fluid-solid equilibrium are in excellent agreement with those reported by Fortini et al.14 in the fluid branch for temperatures up to the vapor-liquid critical point and by Hynninen et al.16 in both branches for temperatures above 0.2. The solid branch is at a lower density than the reported by both works.14,16 See Fig. 5 for a schematic representation of these results. In addition, the vapor-liquid phase diagram obtained by Chapela et al.20 and Fortini et al.14 is also shown. B. Square well potential
Extensive simulations are performed to analyze the vaporliquid and fluid-solid transitions for the SW model with λ = 1.15σ. One phase MC simulations in the NPT ensemble with isotropic fluctuation of volume are performed for nine isobars. Two-phase simulations using the DMD on the NVT ensemble are also carried out to obtain the low density fluidsolid equilibrium using the FCC structure for the solid phase.
FIG. 5. Phase diagram of Yukawa mixtures with κσ = 6. The vapor-liquid coexistence from Chapela et al.20 and Fortini et al.14 is shown with filled squares and open circles, respectively. The fluid-solid coexistence results are also shown: Hynninen et al.16 (filled triangles), Fortini et al.14 (filled circles), and this work using NVT simulations (filled squares).
pressures 0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.5. If the initial conditions of the simulation are within the spinodal region, a violent separation of the liquid and the vapor phases takes place. The vapor-liquid and high density fluid-solid transitions are observed at the lowest pressure. For reduced pressures above or equal 0.5, only the high density-solid transition is observed. Solids with different densities and structures, as in the case of binary Yukawa mixtures, are found. Up to this point, the results look similar to those of long ranged interactions where a low density fluid-solid transition is not observed. The existence of that transition published by several authors12,17,18 is not obtained using NPT simulations. The stability of the solid phase is determined by the pressure tensor components, for example, at T ∗ = 0.58, their values are Py y = 0.226 and Pz z = 0.099. As in the case of the binary colloidal mixtures, they were not the same.
1. SW one phase simulations
The MC results at constant pressure and different temperatures are shown in Fig. 6 for isobars with reduced TABLE II. Components of the pressure tensor for the fluid-solid interfaces obtained from NVT simulations at different temperatures. The total density of the system is also given. T∗ 0.185 0.185 0.185 0.200 0.300 0.400
ρ∗
Px x
Py y
Pz z
0.43 0.68 0.80 0.80 0.80 0.95
0.075 0.122 0.145 0.250 1.020 2.480
0.075 0.122 0.145 0.250 1.020 2.480
0.006 0.006 0.006 0.050 0.800 1.680
FIG. 6. MC NPT results for the SW potential. The fluid-solid and vaporliquid equilibrium envelopes are shown. From bottom to top, the pressures are 0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.5. The filled circles are vapor-liquid curves from this work and open circles correspond to Orkoulas.17 The dashed curve represents the place where the liquid jumps to the solid phase. The low density fluid-solid equilibrium is not observed from NPT simulations.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
054501-6
Méndez-Maldonado et al.
J. Chem. Phys. 142, 054501 (2015)
solids have a CsCl structure, the pressure components are the same. The procedure to obtain the fluid-solid phase diagram should include the analysis of the components of the pressure tensor of the solid phase to establish a mechanical equilibrium necessary condition in addition to the chemical potential, pressure, and temperature. The use of an ideal crystal in the thermodynamic integration method to draw the fluid-solid boundaries leads to a solid with equal diagonal components of the pressure tensor. ACKNOWLEDGMENTS
FIG. 7. Results for the SW using DMD method in the NVT ensemble. Filled squares represent vapor-solid equilibrium curves started with a FCC crystal obtained from runs at densities of 0.641 and 0.9 obtained in this work. Open circles are taken from Orkoulas17 and open triangles from Liu et al.12 The vapor-liquid points obtained in this work are shown with filled circles.
2. SW two-phase simulations
To test the stability of the vapor-liquid phase equilibrium for short ranged interactions, a perturbation from the solution at high density is carried out. If the state is stable, any perturbation must reach the same equilibrium point. The maximum perturbation consists in starting the NVT DMD simulation runs from a solid, obtained from MC NPT simulations. The solid is surrounded by vacuum to obtain a two-phase system. As the system evolves in time, the solid melts and a homogeneous fluid is obtained at the end of the simulation. The reason why the fluid-solid is not formed is because the solid is not stable, as explained above. The crystal structure for the SW model is FCC. When a FCC crystal is surrounded by a low or high density fluid, a low density-crystal equilibrium is obtained. The phase equilibrium results are shown in Fig. 7. The results agree well, within simulation error, with those published by Liu et al.12 and Orkoulas17 using other simulation methods. Small differences at high temperatures are found in the fluid phase where the densities from this work are slightly larger. IV. CONCLUSIONS
Results from this work for short ranged interactions show that the low density fluid-solid equilibrium is only found when the solid has a crystalline structure, CsCl for the binary colloidal mixtures with κσ = 6 and FCC for the square well potential with λ = 1.15σ. The results for binary colloidal mixture using the two-phase simulation method in the NVT ensemble give the same values for the fluid-solid equilibrium as those reported using the free energy calculations14,16 only when the solid phase has a CsCl structure. For the SW fluid, DMD results using the NVT ensemble were around the same as those published by Liu et al.12 and Orkoulas17 using other simulation methods only when the solid has a FCC structure. The solids obtained in NPT simulations at all temperatures have different diagonal components of the pressure tensor, and therefore, they are not mechanically at equilibrium. Therefore, the high density fluid-solid transition is not stable. When the
The authors thank to Laboratorio de Supercómputo for allocation time. J.A. thanks Enrique Velasco for providing us the program to build the CsCl solid structure. We thank the referees for helpful suggestions. 1A.
P. Gast, C. K. Hall, and W. B. Russel, J. Colloid Interface Sci. 96, 251 (1983). 2M. E. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 3991 (1993). 3M. H. J. Hagen and D. Frenkel, J. Chem. Phys. 101(5), 4093 (1994). 4H. N. W. Lekkerkerker, W. C. K. Poon, P. N. Pusey, A. Stroobants, and P. B. Warn, Europhys. Lett. 20, 559 (1992). 5P. N. Pusey, W. C. K. Poon, S. M. Ilett, and P. Barlett, J. Phys.: Condens. Matter 6, A29 (1994). 6E. Leal Calderon, J. Bibette, and J. Biais, Europhys. Lett. 23, 653 (1993). 7E. J. Meijer and D. Frenkel, Phys. Rev. Lett. 67, 1110 (1991). 8G. Foffi, G, D. McCullagh, A. Lawlor, E. Zaccarelli, K. A. Dawson, F. Sciortino, P. Tartaglia, D. Pini, and G. Stell, Phys. Rev. E 65, 031407 (2002). 9L. Mederos and G. Navascués, J. Chem. Phys. 101, 9841 (1994). 10K. P. Shukla, J. Chem. Phys. 112, 10358 (2000). 11M. A. Barroso and A. L. Ferreira, J. Chem. Phys. 116, 7145 (2002). 12H. Liu, S. Garde, and S. Kumara, J. Chem. Phys. 123, 174505 (2005). 13M. Dijkstra, Phys. Rev. E 66, 021402 (2002). 14A. Fortini, A.-P. Hynninen, and M. Dijkstra, J. Chem. Phys. 125, 094502 (2006). 15C. Vega, J. L. F. Abascal, C. McBride, and F. Bresme, J. Chem. Phys. 119, 964 (2003). 16A.-P. Hynninen, M. E. Leunnissen, A. van Blaaderen, and M. Dijkstra, Phys. Rev. Lett. 96, 018303 (2006). 17G. Orkoulas, J. Chem. Phys. 133, 111104 (2010). 18D. L. Pagan and J. D. Gunton, J. Chem. Phys. 122, 184515 (2005). 19J. B. Caballero, A. M. Puertas, A. Fernández-Barbero, F. J. de las Nieves, J. M. Romero-Enrique, and L. F. Rull, J. Chem. Phys. 124, 054909 (2006). 20G. A. Chapela, F. del Río, and J. Alejandre, J. Chem. Phys. 138(5), 054507 (2013). 21M. P. Allen and A. J. Master, Mol. Phys. 79, 277 (1993). 22H. Domínguez, E. Velasco, and J. Alejandre, Mol. Phys. 100, 2739 (2002). 23F. Bresme, M. González-Melchor, and J. Alejandre, J. Phys.: Condens. Matter. 17, S3301 (2005). 24M. González-Melchor, J. Alejandre, and F. Bresme, Phys. Rev. Lett. 90, 135506 (2003). 25M. González-Melchor, F. Bresme, and J. Alejandre, J. Chem. Phys. 122, 104710 (2005). 26J. Alejandre, F. Bresme, and M. González-Melchor, Mol. Phys. 107, 357 (2009). 27M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, J. Phys. A: Math. Gen. 39, 5629 (2006). 28M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulations (Oxford University Press, New York, 2010). 29D. Frenkel and B. Smith, Understanding Molecular Simulation from Algorithms to Applications, 2nd ed. (Academic Press, San Diego, 1996). 30P. Gregory and J. C. Schug, Mol. Phys. 82(4), 677 (1994). 31B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, 459 (1959). 32D. C. Rapaport, J. Chem. Phys. 71, 3299 (1979). 33G. A. Chapela, S. E. Martínez-Casas, and J. Alejandre, Mol. Phys. 53, 139 (1984). 34P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford University Press, New York, 1987).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 15:38:28
