Phase diagram of colloidal systems: A comparison of standard monte carlo simulations and virtual-Move monte carlo simulations

HNUE JOURNAL OF SCIENCE DOI: 10.18173/2354-1059.2018-0032 Natural Science, 2018, Volume 63, Issue 6, pp. 84-91 This paper is available online at PHASE DIAGRAM OF COLLOIDAL SYSTEMS: A COMPARISON OF STANDARD MONTE CARLO SIMULATIONS AND VIRTUAL-MOVE MONTE CARLO SIMULATIONS Pham Van Hai Faculty of Physics, Hanoi National University of Education Abstract. In this work, we investigate the phase diagram of a three-dimensional Lennard-Jones (LJ) and square-well (SQW) system by means of both standardMonte Carlo simulation (SPMC) and “virtual-move” Monte Carlo simulations (VMMC). We find that the phase diagrams obtained by both SPMC and VMMC are very similar. However, there exist a few significant differences in the phase formation process as well as in the final phase states, especially in cases of large interaction energies. We interpret the discrepancy between the two methods as results of a slow equilibrium process in SPMC, in contrast to a more rapid equilibration in VMMC. These findings, in good agreement with those observed through the self-assembly of strongly attractive systems in the literature, suggests that the VMMC method produces more accurate phase diagrams. Keywords: Lennard-Jones potential, square-well potential, virtual move Monte Carlo, phase diagram. 1. Introduction Colloidal suspension is a mixture of a dispersed solid with particle sizes from micrometers to nanometers and a continuous liquid phase. Colloidal suspension, such as milk, gelatin and wet clay, plays important roles in pharmaceuticals, food, non-drip paints, sewage disposal, and wet photography [1]. The phase diagram of colloidal systems can be well-described by a range of simple model systems, including a hard-sphere model [2, 3], repulsive screened Coulomb (Yukawa) model [4, 5], Lennard-Jones model [6, 7, 8] and square-well model [9, 10]. However, in all of these standard Monte Carlo methods employed, each simulation run only consists of a subsequent movement of single individual particles and neglect collective modes, i.e. rotational and translational movement of clusters are absent [11]. This could lead to kinetic traps as result of a low acceptance rate, especially in colloidal systems with a strong attractive inter-particle interaction [12, 13]. To avoid unphysical kinetics traps in MC simulations of systems with highly attractive interactions, Whitelam and Geissler [12, 13] have introduced a modified MC algorithm Received July 25, 2018. Revised August 8, 2018. Accepted August 15, 2018. Contact Pham Van Hai, e-mail: haipv@hnue.edu.vn. 84 Phase diagram of colloidal systems: A comparison of standard monte carlo simulations and virtual-move... where virtual moves of clusters are taken into account. However, the authors have only considered a simple case of two-dimensional Lennard-Jones disk for several points in the phase diagram. Our aim is here to employ the MC algorithm by Whitelam and Geissler [12] to calculate the full phase diagrams for some familiar systems. At the same time, we compare the obtained results with those derived from standard Monte Carlo simulation. 2. Content 2.1. Model and methods Figure 1. Simple isotropic potential models for colloidal systems (a) Lennard-Jones potential and (b) square-well potential Consider a simple model of a system for which the pair potential U(r) at a separation r gives ULJ(r) = 4ǫ [( r σ )12 − ( r σ )6] , (2.1) Equation (2.1) is the famous 12-6 potential of Lennard-Jones that involves two parameters: ǫ, the depth of the potential well the minimum value in U(r)); and σ, the collision diameter at which the separation of the particles U(r) = 0. The Lennard-Jones potential is widely used to describe the interaction between rare-gas atoms and spherical-like molecules [14, 15], as illustrated in Figure 1a. The square-well potential is also determined by two parameters, ǫ and λ; ǫ is the depth of the well, and (λ− 1)σ is the width of the well. The potential shape with a typical value λ of 1.2 is depicted in Figure 1b. Despite its simplicity, the model of the square-well potential has some applications to real, physical systems, e.g. a suspension of micron-sized silica spheres dispersed in an organic solvent [15]. The explicit form of the square-well potential is following: USQW(r) =   ∞ r < σ −ǫ σ < r < λσ 0 otherwise (2.2) 85 Pham Van Hai We carry out all simulations in the canonical ensemble (NV T ensemble) using both the SPMC and VMMC method. In the SPMC method [16, 17, 18], each simulation run includes initialization, equilibration and sampling. After initialization, the system is evolved with MC steps. Each MC step consists of N (N is the number of particles) attempts to randomly displace a single particle. This trial move either accepted or rejected according to a Metropolis probability of 50%. The virtual-move Monte Carlo simulation also uses the Metropolis algorithm. However, this method permits of both realistic motion and collective internal rearrangement. The collective movement drives the formation of kinetics traps and facilitates orderly growth. The good approximations of a realistic dynamics apply for strong, short-ranged interaction and a high degree of angular specificity. In a short interaction range, the affection of interaction energy of a particle on others is higher than in long-ranged of interaction. The algorithm of VMMC method is given in detail in Ref [12]. In general, a virtual trial move of a particle is proposed. Then the bond energy of the particle and its neighbors is calculated before and after the trial moves are accepted. The procedure for all neighbors is iterated until no more particles moves. The positions of all affected particles are updated simultaneously. In our simulations, we perform with N = 500 particles that are initialized randomly in the cubic box with the periodic boundary condition. The total number of MC cycles is 106. For a given set of parameters, five independent simulation runs are performed to improve statistical accuracy. 2.2. Result and discussion Figure 2. (Left panel) Number density as a function of the scaled distance for two different types of phases. For the sake of clarity, here we assume that density profile is distributed along the z-axis (Right panel) Radial distribution function for solid phase and liquid phase We determine the phase diagrams in the plane η − T ∗, where η is the particle density and T ∗ is the reduced temperature defined as T ∗ = kBT/ǫ. In the Lennard-Jones potential, we choose the number density range from 0.1 to 0.7 in the step of 0.2 and the reduced temperature in range of 0.2-2.0 in the step of 0.3. In case of square-well potential, we choose the number density ranging from 0.2 − 0.7 in the step of 0.1 and temperature in the range of 0.1 − 1.3 with the step of 0.1. In order to determine whether the state points in the phase diagram are gas, liquid, solid, or phase coexistence, we employ the density profile, radial distribution function, in combination with 86 Phase diagram of colloidal systems: A comparison of standard monte carlo simulations and virtual-move... visual inspection. First, we use the visual snapshots of systems under consideration [20] to guess whether there exists a single phase or the phase separation. If the system stays in a gas phase, the movement of particles is fully chaotic. Furthermore, the motion of particles is much faster than that of particles in the liquid phase or solid phase. If the particle belongs to the liquid phase, it moves around but keeping connections with other particles. In other words, a particle moves by sliding with others. If the particle belongs to the solid phase, we observe an order arrangement of particles. To distinguish the phase separation, including gas-liquid, gas-solid, and liquid-solid, one can use the density profile [19]. The density profile of a homogeneous system fluctuates around its number density, whereas the histogram of a system that phase separates has two maxima at the number densities of the two separated areas, respectively (see Figure 2a for illustration). Hence this enables us not only to detect phase separation, but to read off the packing fractions of the coexisting phases as well. A very important means to analyze the structure of a system is via its RDF g(r), which is defined as the fraction of the average number density at a distance r from a given particle and the number density of an ideal gas with the same overall density. This function can be used to calculate further important observables of the system such as the average potential energy and the pressure. Furthermore, g(r) is important because it contains information about the local structure of the system. In particular fluid and solid phases can be distinguished by means of the RDF, as shown in Figure 2b. 2.2.1. Qualitative three-dimensional phase diagram of Lennard-Jones potential Figure 3. Phase diagrams of a three-dimensional Lennard-Jones (LJ) system in the reduced temperature T ∗-number density η representation using (a) SPMC and (b) VMMC method In Figure 3, we show comparative phase diagrams obtained by SPMC and VMMC, respectively. It can be seen that two diagrams are very similar to each other, in agreement with those reported previously [6, 8]. They made of gas phase, liquid and gas-liquid coexistence at T ∗ > 0.6 and η = 0.1 − 0.7. However, different from Ref. [14] where the authors observed gas-solid coexistence at T ∗ = 0.5, we find a gas-liquid coexistence. Furthermore, at T ∗ = 0.2, 0.5, and 0.7, a gas-solid coexistence is observed in simulations using the SPMC method (marked by red circle) instead of a homogenous liquid phase using the VMMCmethod, as shown by visual 87 Pham Van Hai snapshots in Figure 4. In addition, the VMMC method enables the equilibration rate much more rapid than that of the SPMC, for example T ∗ = 0.2, η = 0.5 the SPMC requires 106 MC cycles, while the VMMC needs only 6.35 × 105 MC cycles to reach an equilibrium state. Figure 4. Simulation snapshot obtained at T ∗ = 0.2, η = 0.5, by (a) SPMC method and (b) VMMC method 2.2.2. Qualitative three-dimensional phase diagram of square-well potential Figure 5. Phase diagrams of a three-dimensional square-well potential (SQW) system using SPMC (a) and VMMC (b) We performed computer simulations for the SQW potential for three cases of the interaction range of λ, but only the λ = 2 case is shown here as an example. Figure 5 shows the representative phase diagram in the T ∗− η plane for the SQW potential. We note that it is difficult to distinguish the phase separation into a gas-liquid region and gas-solid region in the SQW model due to 88 Phase diagram of colloidal systems: A comparison of standard monte carlo simulations and virtual-move... a disordered arrangement of spheres in the solid phase. The imperfect phase separation could be gas-solid coexistence with meta-stable fluid phase. Therefore, in Figure 5 the state points corresponding to metastable phases are also shown. The difference in the phase diagram between the SPMC and VMMC method lies mostly at T ∗ = 1.2, η = 0.7 and T ∗ = 1.2, η = 0.5. For clarity, we shows representative simulation snapshots based on the SPMC method and VMMC method at these two points (see Figure 6). In addition, the phase formation process occurs differently between two methods. The equilibrium state obtained by SPMC takes a longer MC time compared to that of VMMC. Figure 7 illustrates the visual snapshots taken at different stages of the simulation process. Apparently, at T ∗ = 1.2, η = 0.7 SPMC needs 5.24 × 105 MC steps but the VMMC just requires 2 × 104 MC steps for equilibration. The explanation as follows: the VMMC method takes into account the collective motion mode of clusters, the system seems to move rapidly more than simulation by using the SPMC method. Figure 6. Simulation snapshots at T ∗ = 1.2, η = 0.7 obtained by SPMC (a) and VMMC (b), at T ∗ = 1.2, η = 0.5 obtained by SPMC (c) and VMMC (d) 89 Pham Van Hai Figure 7. Simulation snapshots taken at T ∗ = 1.2, η = 0.7 in case of figure (a1) initial configuration, (a2) state after 5.24 × 105 MC steps and (a3) final state of the system simulated by SPMC method, while figure (b1) represents the initial state, figure (b2) is state after 2× 104 MC steps and figure (b3) is the final state of the system using VMMC method 3. Conclusion In summary, we have investigated the phase diagram for one-component systems with Lennard-Jones and Square-Well interacting potential. In each case, we performed the standard MC simulations and “virtual-move” MC simulations for comparative purpose. Some main conclusions can be draw as follows. In almost case of the points at η = 0.1 − 0.7 in the phase diagram, our results show good agreement with those obtained by previous experiments and other computer simulations. The interesting discrepancy is at reduced temperature T ∗ = 0.5. For the case of the Lennard-Jones potential, the SPMCmethod produces the gas-solid coexistence, while the VMMC method favors the homogenous liquid phase at T ∗ = 0.5 and η in the range from 0.5 to 0.7. Similar results achieved in the case of the square-well potential in which we find some difference in the resulting phase state. For example, a homogenous crystal phase is presented in the system using SPMC instead of a liquid phase in the system using VMMC. In particular, a discrepancy between two simulation methods increases with increasing of the interaction energy. In addition, from visual snapshots we observed that the equilibrium state obtained by SPMC takes a longer MC time compared to that of VMMC. We interpret this as the collective motion modes taken into account the VMMC method. This also suggests that in the same intervals, the VMMC method could produce the equilibrium state more accurate than the method of SPMC method. Acknowledgment. This work was funded byVietnamese National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.02-2017.328. 90 Phase diagram of colloidal systems: A comparison of standard monte carlo simulations and virtual-move... REFERENCES [1] T. Cosgrove, 2005. Colloid Science: Principles, Methods and Applications. Blackwell Publishing. [2] B. J. Alder and T. E. Wainwright, 1957. Journal of Chemical Physics. 27, 1208. [3] Z. Cheng, P.M. Chaikin, W.B. Russel, W.V. Meyer, J. Zhu, R.B. Rogers, and R.H. Ottewill, 2001. Materials & Design. 22, 529. [4] Kurt Kremer, Mark O. Robbins, and Gary S. Grest, 2001. Physical Review Letter. 57, 2694. [5] M. O. Robbins, 1988. Journal of Chemical Physics. 88, 3286. [6] B. Smit, 1992. Journal of Chemical Physics. 96, 8639. [7] S.T Lin, M. Blanco, and W.A. Goddard, 2003. Journal of Chemical Physics. 119, 11792. [8] H. Watanabe, N. Ito, and C.K. Hu, 2012. Journal of Chemical Physics. 136, 204102. [9] D. L. Pagana and J. D. Gunton, 2005. Journal of Chemical Physics. 122, 184515. [10] J.C. Armas-Pérez, J. Quintana, G.A. Chapela, E. Velasco, and G. Navascués, 2014 Journal of Chemical Physics 140, 064503. [11] S. Ru˚zˇicˇka and Michael P. Allen 2014. Physical Review E 89, 033307. [12] S. Whitelam and P.L Geissler. 2007. Journal of Chemical Physics 127, 154101. [13] S. Whitelam. 2011. Molecular Simulation 37, 606. [14] J.P. Hansen and L. Verlet, 1969. Physical Review 184, 151. [15] Jean-Pierre Hansen, I.R. McDonald, 2013. Theory of Simple Liquids: with Applications to Soft Matter. Academic Press. [16] Hai P.V., A. Fortini and M. Schmidt 2016. Physical Review E 93, 052609. [17] Hai P.V., A. Fortini and M. Schmidt 2017. Materials 10, 11. [18] Hai P.V., Hoang L.H. and Thuy N.M. 2018. Colloid and Polymer Science 296, 1387. [19] T. Giorgino. 2014. Computer Physics Communications 185, 317. [20] W. Humphrey, A. Dalke, and K. Schulten, 1996. Journal of Molecular Graphics 14, 33. 91