(168b) Precise Simulation of Fluid-Solid Transitions Using Constrained Cell Models

In simulation studies of fluid-solid transitions, the constrained cell model is defined by dividing the simulation volume, V, into N Wigner-Seitz cells appropriate for the solid phase under consideration. Each of the N particles is assigned and confined to move in a single Wigner-Seitz cell of volume V/N. The constrained cell model is a realistic representation of the high-density solid phase. In the analysis of the freezing transition of a system of hard spheres via thermodynamic integration, Hoover and Ree [1, 2] used the constrained cell model to obtain the free energy of the solid phase. The simulations indicated that the pressure vs. density curve of the solid (as modeled through the constrained cell model) exhibits a cusp (or kink or angular point) at a density which is approximately 64% of the density at close packing. The appearance of such an anomaly was related to the mechanical stability of the solid. Specifically, for densities that are lower than the density associated with the cusp, the solid phase cannot survive without the presence of the walls of the Wigner-Seitz cells and it quickly transforms to a disordered, fluid-like phase. At higher densities, the solid phase can exist for substantial time intervals without the presence of the cell walls. Although the presence of such an instability or anomaly does not affect free energy calculations, it is the main reason for which the constrained cell model has not been very popular in simulation studies of fluid-solid transitions.

In thermodynamic integration techniques based on constrained cell models [2], one usually extends the simulations from the high-density, nearly-incompressible limit to the low-density, ideal-gas region for which the free energy can be obtained analytically. In order to reduce the number of simulated states, Hoover and Ree proposed [1] a more general cell model by adding a homogenous external field that controls the relative stability of the solid versus the fluid phase. High field values force single occupancy configurations with one particle per Wigner-Seitz cell and thus favor the solid phase. Normal (unconstrained) system behavior is restored in the limit of vanishing field. Hence, the constrained cell model is a special case of the more general or modified cell model. Hoover and Ree thought that the modified cell model could be used to link the fluid with the solid phase on a constant-density (or pressure) path by gradually increasing the strength of the field, thus reducing the number of simulated states in thermodynamic integration techniques. In their original work associated with hard-sphere freezing, Hoover and Ree investigated the properties of this model at low densities through cluster expansion techniques [1].  This model has received very little (if any) attention in simulations of freezing transitions.

In the present work, the modified cell model is reformulated in the isothermal-isobaric ensemble and its use in simulation studies of fluid-solid transitions is considered [3-7]. The simulation results indicate that as the strength of the field is reduced at constant pressure, the transition from the solid to the fluid phase is smooth and continuous at low and moderate pressures. In contrast, at high pressures, the passage from the solid to the fluid occurs via a discontinuous first-order phase transition. The special point that separates continuous from discontinuous behavior is very close to the mechanical stability point, which now appears as an inflection point in the pressure-density isotherm of the solid phase. Hence, the mechanical stability point plays the role of some type of critical/Curie point for the modified cell model. The fluid-solid transition of the Lennerd-Jones model system on an isotherm that corresponds to a reduced temperature of 2 is determined by analyzing the field-induced, order-disorder transition of the corresponding modified cell model in the high-pressure, vanishing-field limit using flat-histogram techniques [5]. It is also shown that this model can be used to link the fluid with the solid phase, thus reducing the number of simulated states in thermodynamic integration [6].

The use of the modified cell model in simulations of fluid-solid transitions on subcritical isotherms is also considered [7]. Specifically, the fluid-solid transition of a system of Lennard-Jones particles is determined for the isotherm that corresponds to a reduced temperature of unity, which is lower than the known critical temperature of this system. Constant-pressure simulations of the corresponding constrained cell model at a reduced temperature of unity show that there is a first-order phase transition between a dilute, gas-like and a dense, liquid-like phase, respectively. This transition, which is the analogue of the gas-liquid transition of the unconstrained system, terminates at a critical point. The corresponding phase diagram of the constrained cell model is determined through constant-pressure simulations. It is found that the critical temperature and pressure of the constrained system are both higher than those of the unconstrained system. The main reason for this behavior is the reduction of the entropy of the system caused by the single occupancy constraint. Apart from this first-order, gas-liquid phase transition, the behavior of the pressure-density subcritical isotherms of the fluid (unconstrained) and the solid (constrained) phase is similar to supercritical isotherms. The fluid-solid transition of the Lennard-Jones model at a reduced temperature of unity is obtained by analyzing the field-induced phase transition of the modified cell model and by thermodynamic integration using the same model [7].

In all cases considered in this work [3-7], fluid-solid coexistence is determined through finite-size scaling techniques for first-order phase transitions. These scaling techniques are based on analyzing the size-dependent behavior of susceptibilities and dimensionless moment ratios of the order parameter. The data for small systems appear to be irregular and they do not conform at all to the expected scaling forms [5-7]. The main reason for this peculiar behavior is the approximate nature of accounting for long-range corrections in the energy evaluation, which is more pronounced for small systems. Since these effects become progressively smaller as the system size increases, the results for moderate and large systems conform well to the scaling laws and yield reliable estimates of phase coexistence in the limit of infinite size. The scaling analysis clearly demonstrates the importance of accounting for size effects in simulation studies of fluid-solid transitions.

[1] W. G. Hoover and F. H. Ree, J. Chem. Phys. 47, 4873 (1967).

[2] W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968).

[3] G. Orkoulas and M. Nayhouse, J. Chem. Phys. 134, 171104 (2011).

[4] M. Nayhouse, A. M. Amlani, and G. Orkoulas, J Phys.: Condens. Matter 23, 325106 (2011).

[5] M. Nayhouse, A. M. Amlani, and G. Orkoulas, J. Chem. Phys. 135, 154103 (2011).

[6] M. Nayhouse, A. M. Amlani, V. R. Heng, and G. Orkoulas, J. Phys.: Condens. Matter 24, 155101 (2012).

[7] M. Nayhouse, V. R. Heng, A. M. Amlani, and G. Orkoulas, J. Phys. A: Mathematical and Theoretical 45, 155002 (2012).