Solid-Fluid Equilibrium: Insights from Simple Molecular Models. In: Prigogine I, Rice SA, editors. Advances in Chemical Physics. Hoboken: John Wiley & Sons, Inc.; 2000. pp. 113-79. [DOI: 10.1002/9780470141748.ch2]

Efficient approaches to solutions of partition function for condensed matters

The key problem of statistical physics standing over one hundred years is how to exactly calculate the partition function (or free energy), which severely hinders the theory to be applied to predict the thermodynamic properties of condensed matters. Very recently, we developed a direct integral approach (DIA) to the solutions and achieved ultrahigh computational efficiency and precision. In the present work, the background and the limitations of DIA were examined in details, and another method with the same efficiency was established to overcome the shortage of DIA for condensed system with lower density. The two methods were demonstrated with empirical potentials for solid and liquid cooper, solid argon and C₆₀ molecules by comparing the derived internal energy or pressure with the results of vast molecular dynamics simulations, showing that the precision is about ten times higher than previous methods in a temperature range up to melting point. The ultrahigh efficiency enables the two methods to be performed with ab initio calculations and the experimental equation of state of solid copper up to ∼600 GPa was well reproduced, for the first time, from the partition function via density functional theory implemented.

Monte Carlo simulation and integral equation study of Hertzian spheres in the low-temperature regime

We investigate the behavior of Hertzian spheres in the fluid phase and in proximity of the freezing threshold by using Monte Carlo (MC) simulations and integral equation theories, based on the Ornstein-Zernike (OZ) approach. The study is motivated by the importance of the Hertzian model in representing a large class of systems interacting via soft interactions, such as star polymers or microgels. Radial distribution functions, structure factors, and excess entropy clearly show the reentrant behavior typical of the Hertzian fluid, well caught by both MC simulations and OZ theory. Then, we make use of some phenomenological one-phase criteria for testing their reliability in detecting the freezing threshold. All criteria provide evidence of the fluid-solid transition with different degrees of accuracy. This suggests the possibility to adopt these empirical rules to provide a quick and reasonable estimate of the freezing transition in model potentials of wide interest for soft matter systems.

Study on the Thermodynamics of Polymer Crystallization Based on Twin-Lattice Model

Polymer crystallization is the most important part in determining the performance of polymeric materials. The twin-lattice model originally provided by Lennard-Jones and Devonshire, developed by Pople and Karasz and other researchers, is extended for describing the thermodynamics of polymer crystallization. The positional order of segments and the orientational order of bonds are considered in this model. The free energy of polymers is obtained by further introducing the conformational energy and entropy, and thus a new parameter is defined, which is the ratio of conformational energy and positional diffusion energy. We studied two kinds of processes in polymer crystallization, including the process with plastic crystal phase and without any mesophases. The choice of crystallizing process is determined by the magnitude of lattice energy and conformational energy. The solid-solid transition from crystal to plastic crystal shows a significant dependence on the conformational energy. Considering data reliability, n-paraffins are chosen as the representation of polymers to compare the predictions of the model with experimental observations. We predict the number of carbons beyond which the rotator phase disappears, which is quite in agreement with the experiments. These calculations and results show this model can provide a new understanding to the crystallization of polymers.

Efficient approach to the free energy of crystals via Monte Carlo simulations

We present a general approach to compute the absolute free energy of a system of particles with constrained center of mass based on the Monte Carlo thermodynamic coupling integral method. The version of the Frenkel-Ladd approach [J. Chem. Phys. 81, 3188 (1984)]JCPSA60021-960610.1063/1.448024, which uses a harmonic coupling potential, is recovered. Also, we propose a different choice, based on one-particle square-well coupling potentials, which is much simpler, more accurate, and free from some of the difficulties of the Frenkel-Ladd method. We apply our approach to hard spheres and compare with the standard harmonic method.

Crystal growth kinetics in Lennard-Jones and Weeks-Chandler-Andersen systems along the solid-liquid coexistence line

Kinetics of crystal-growth is investigated along the solid-liquid coexistence line for the (100), (110), and (111) orientations of the Lennard-Jones (LJ) and Weeks-Chandler-Andersen (WCA) fcc crystal-liquid interface, using non-equilibrium molecular dynamics simulations. A slowing down of the growth kinetics along the coexistence line is observed, which is due to the decrease of the melting enthalpy with increasing coexistence temperature and pressure. Other quantities such as the melting pressure and liquid self-diffusion coefficient have a comparatively lesser impact on the kinetic growth coefficient. Growth kinetics of the LJ and WCA potentials become similar at large values of the melting temperature and pressure, when both resemble a purely repulsive soft-sphere potential. Classical models of crystallization from the melt are in reasonable qualitative agreement with our simulation data. Finally, several one-phase empirical melting/freezing rules are studied with respect to their validity along the coexistence line.

Effect of confinement on the solid-liquid coexistence of Lennard-Jones fluid

The solid-liquid coexistence of a Lennard-Jones fluid confined in slit pores of variable pore size, H, is studied using molecular dynamics simulations. Three-stage pseudo-supercritical transformation path of Grochola [J. Chem. Phys. 120(5), 2122 (2004)] and multiple histogram reweighting are employed for the confined system, for various pore sizes ranging from 20 to 5 molecular diameters, to compute the solid-liquid coexistence. The Gibbs free energy difference is evaluated using thermodynamic integration method by connecting solid-liquid phases under confinement via one or more intermediate states without any first order phase transition among them. Thermodynamic melting temperature is found to oscillate with wall separation, which is in agreement with the behavior seen for kinetic melting temperature evaluated in an earlier study. However, thermodynamic melting temperature for almost all wall separations is higher than the bulk case, which is contrary to the behavior seen for the kinetic melting temperature. The oscillation founds to decay at around H = 12, and beyond that pore size dependency of the shift in melting point is well represented by the Gibbs-Thompson equation.

Algorithm for constant-pressure Monte Carlo simulation of crystalline solids

We describe an alternative method for performing isothermal-isobaric Monte Carlo simulations of crystalline solids. This method uses thermodynamics to estimate appropriate scaling of coordinates relative to their nominal lattice sites in order to increase the probability of acceptance of volume changes. We test this coordinate scaling with three systems: hard spheres, Lennard-Jones spheres, and hard dumbbells. In all cases, we find that the move allows for both a larger step size and faster convergence of calculated properties, in comparison to the conventional algorithm. The improvement is more dramatic for hard potentials, where compressing the system using the conventional algorithm (even a small amount) will almost always cause an overlap.

Universal features of the free-energy functional at the freezing transition for repulsive potentials

The free-energy difference between coexisting solid and liquid phases is studied in the context of classical density functional theory (DFT). A bridge function is used to represent the higher-order (n>2) terms in the perturbative expansion of the excess Helmholtz free energy, and the values of this bridge function within the solid lattice are determined by inversion using literature Monte Carlo simulation results. Four potential models, specifically hard-sphere and inverse twelfth-, sixth-, and fourth-power repulsive, are studied. The face-centered cubic (fcc) solid is considered for the hard-sphere and inverse twelfth- and sixth-power potentials, while the body-centered cubic (bcc) solid is considered for the inverse sixth- and fourth-power potentials. For a given solid structure there is a remarkable similarity among the bridge functions for different potentials that is analogous to the universality in the sum of elementary diagrams, or bridge functions, of liquid-state theory as originally observed by Rosenfeld and Ashcroft [Phys. Rev. A 20, 1208 (1979)]. In further analogy with liquid-state theory, the bridge functions in the present problem are plotted as functionals of the second-order convolution term in the perturbative expansion. In each case, the plot indicates a unique functionality in the dense regions of the solid near the lattice sites but a scattered and nonunique behavior in the void regions. Interestingly, knowledge of the functional relationship in the unique region near the lattice sites seems to be sufficient to quantitatively model the solid-fluid phase transition. These qualitative observations are true for both fcc and bcc solid phases, although there are some quantitative differences between them. The findings suggest that pursuit of a closure-based DFT of solid-fluid transitions may be profitable.

Excess entropy and structural transitions in a two-dimensional square-shoulder fluid

Metropolis Monte Carlo simulations on the square-shoulder fluid of Malescio and Pellicane are used to trace the temperature dependent excess entropy, the heat capacity, and configurational energy along several isochores, including those for which mechanically stable zero-temperature structures have been predicted. Thermodynamic signatures of structural phase transitions are identified along several isochores, in addition to the low-density triangular solid and stripe phase transitions identified earlier. The finite temperature phases illustrate the competition between cluster formation and stripe formation as competing mechanisms for generating minimum free energy configurations as a function of density, consistent with earlier results at zero temperature. We also critically examine the usefulness of a phase-ordering rule based on the residual multiparticle entropy (RMPE) in predicting the formation of this diverse set of ordered structures from a disordered fluid phase. For the majority of the isochores studied, the RMPE prediction and the thermodynamic evidence for a phase transition were consistent. However, this criterion fails along isochores that are in regions of coexistence. Thus, the zero-RMPE rule is only likely to be approximately predictive in systems with small phase coexistence regimes, e.g., in the case of liquid crystal forming systems.

Residual multiparticle entropy does not generally change sign near freezing

The residual multiparticle entropy (RMPE) of two- and three-dimensional fluids changes sign near the freezing line, providing a quasiuniversal "one-phase" rule for the location of the liquid-solid transition. We present new simulation results for d-dimensional hard-sphere fluids (d=1-5) which show, however, that this freezing criterion fails in other spatial dimensions. The results also call into question the idea that a change in sign of the RMPE implies the emergence of a new kind of local structural order in the fluid.

Vapor-liquid equilibria from the triple point up to the critical point for the new generation of TIP4P-like models: TIP4P/Ew, TIP4P/2005, and TIP4P/ice

The vapor-liquid equilibria of three recently proposed water models have been computed using Gibbs-Duhem simulations. These models are TIP4P/Ew, TIP4P/2005, and TIP4P/ice and can be considered as modified versions of the TIP4P model. By design TIP4P reproduces the vaporization enthalpy of water at room temperature, whereas TIP4P/Ew and TIP4P/2005 match the temperature of maximum density and TIP4P/ice the melting temperature of water. Recently, the melting point for each of these models has been computed, making it possible for the first time to compute the complete vapor-liquid equilibria curve from the triple point to the critical point. From the coexistence results at high temperature, it is possible to estimate the critical properties of these models. None of them is capable of reproducing accurately the critical pressure or the vapor pressures and densities. Additionally, in the cases of TIP4P and TIP4P/ice the critical temperatures are too low and too high, respectively, compared to the experimental value. However, models accounting for the density maximum of water, such as TIP4P/Ew and TIP4P/2005 provide a better estimate of the critical temperature. In particular, TIP4P/2005 provides a critical temperature just 7 K below the experimental result as well as an extraordinarily good description of the liquid densities from the triple point to the critical point. All TIP4P-like models present a ratio of the triple point temperature to the critical point temperature of about 0.39, compared with the experimental value of 0.42. As is the case for any effective potential neglecting many body forces, TIP4P/2005 fails in describing simultaneously the vapor and the liquid phases of water. However, it can be considered as one of the best effective potentials of water for describing condensed phases, both liquid and solid. In fact, it provides a completely coherent view of the phase diagram of water including fluid-solid, solid-solid, and vapor-liquid equilibria.

Entropy, local order, and the freezing transition in Morse liquids

The behavior of the excess entropy of Morse and Lennard-Jones liquids is examined as a function of temperature, density, and the structural order metrics. The dominant pair correlation contribution to the excess entropy is estimated from simulation data for the radial distribution function. The pair correlation entropy (S2) of these simple liquids is shown to have a threshold value of (-3.5+/-0.3)kB at freezing. Moreover, S2 shows a T(-2/5) temperature dependence. The temperature dependence of the pair correlation entropy as well as the behavior at freezing closely correspond to earlier predictions, based on density functional theory, for the excess entropy of repulsive inverse power and Yukawa potentials [Rosenfeld, Phys. Rev. E 62, 7524 (2000)]. The correlation between the pair correlation entropy and the local translational and bond orientational order parameters is examined, and, in the case of the bond orientational order, is shown to be sensitive to the definition of the nearest neighbors. The order map between translational and bond orientational order for Morse liquids and solids shows a very similar pattern to that seen in Lennard-Jones-type systems.

Determining landscape-based criteria for freezing of liquids

The correlation between statistical properties of the energy landscape and the number of accessible configurational states, as measured by the exponential of the excess entropy (e(Se)), are studied in the case of a simple Lennard-Jones-type liquid in the neighborhood of the thermodynamic freezing transition. The excess entropy Se is defined as the difference between the entropy of the liquid and that of the ideal gas under identical temperature and pressure conditions and is estimated using the pair correlation contribution, S2. Landscape properties associated with three categories of configurations are considered: instantaneous configurations, inherent saddles, and inherent minima. Landscape properties studied include the energy and the key parameters of the Hessian eigenvalue distribution as well as the mean distances between instantaneous configurations and the corresponding inherent saddles and minima. The signatures of the thermodynamic freezing transition are clearest in the case of inherent structure properties which show, as a function of e(S2), a pronounced change in slope in the vicinity of the solid-liquid coexistence. The mean distance between instantaneous and saddle configurations also shows a similar change in slope when the system crosses from the stable to the supercooled regime. In the case of inherent saddles, the minimum eigenvalue acts as a similar indicator of the thermodynamic freezing transition but the average and maximum eigenvalues do not carry similar signatures. In the case of instantaneous configurations, a weak indicator of the thermodynamic freezing transition is seen in the behavior of the fraction of negative curvature directions as a function of the exponential of the excess entropy.

Further studies of a simple atomistic model of silica: Thermodynamic stability of zeolite frameworks as silica polymorphs

We have applied our previously reported model of silica based on low coordination and strong association [J. Chem. Phys. 121, 8415 (2004)], to the calculation of phase stability of zeolite frameworks SOD, LTA, MFI, and FAU as silica polymorphs. We applied the method of Frenkel and Ladd for calculating free energies of these solids. Our model predicts that the MFI framework structure has a regime of thermodynamic stability at low pressures and above approximately 1400 K, relative to dense phases such as quartz. In contrast, our calculations predict that the less dense frameworks SOD, LTA, and FAU exhibit no regime of thermodynamic stability. We have also used our model to investigate whether templating extends the MFI regime of thermodynamic stability to lower temperatures, by considering templates with hard-sphere repulsions and mean-field attractions to silica. Within the assumptions of our model, we find that quartz remains the thermodynamically stable polymorph at zeolite synthesis temperatures (approximately 400 K) unless unphysically large template-silica attractions are assumed. These predictions suggest that some zeolites such as MFI may have regimes of thermodynamic stability even without template stabilization.

Crystal nucleation in binary hard sphere mixtures: A Monte Carlo simulation study

We present calculations of the nucleation barrier during crystallization in binary hard sphere mixtures under moderate degrees of supercooling using Monte Carlo simulations in the isothermal-isobaric semigrand ensemble in conjunction with an umbrella sampling technique. We study both additive and negatively nonadditive binary hard sphere systems. The solid-fluid phase diagrams of such systems show a rich variety of behavior, ranging from simple spindle shapes to the appearance of azeotropes and eutectics to the appearance of substitutionally ordered solid phase compounds. We investigate the effect of these types of phase behavior upon the nucleation barrier and the structure of the critical nucleus. We find that the underlying phase diagram has a significant effect on the mechanism of crystal nucleation. Our calculations indicate that fractionation of the species upon crystallization increases the difficulty of crystallization of fluid mixtures and in the absence of fractionation (azeotropic conditions) the nucleation barrier is comparable to pure fluids. We also calculate the barrier to nucleation of a substitutionally ordered compound solid. In such systems, which also show solid-solid phase separation, we find that the phase that nucleates is the one whose equilibrium composition is closer to the composition of the fluid phase.

Evaluation of phenomenological one-phase criteria for the melting and freezing of softly repulsive particles

We test the validity of some widely used phenomenological criteria for the localization of the fluid-solid transition thresholds against the phase diagrams of particles interacting through the exp-6, inverse-power-law, and Gaussian potentials. We find that one-phase rules give, on the whole, reliable estimates of freezing/melting points. The agreement is ordinarily better for a face-centered-cubic solid than for a body-centered-cubic crystal, even more so in the presence of a pressure-driven reentrant transition of the solid into a denser fluid phase, as found in the Gaussian-core model.

Freezing line of the Lennard-Jones fluid: A phase switch Monte Carlo study

We report a phase switch Monte Carlo (PSMC) method study of the freezing line of the Lennard-Jones (LJ) fluid. Our work generalizes to soft potentials the original application of the method to hard-sphere freezing and builds on a previous PSMC study of the LJ system by Errington [J. Chem. Phys. 120, 3130 (2004)]. The latter work is extended by tracing a large section of the Lennard-Jones freezing curve, the results of which we compare with a previous Gibbs-Duhem integration study. Additionally, we provide new background material regarding the statistical-mechanical basis of the PSMC method and extensive implementation details.

Thermodynamic properties of lattice hard-sphere models

Thermodynamic properties of several lattice hard-sphere models were obtained from grand canonical histogram- reweighting Monte Carlo simulations. Sphere centers occupy positions on a simple cubic lattice of unit spacing and exclude neighboring sites up to a distance sigma. The nearestneighbor exclusion model, sigma = radical2, was previously found to have a second-order transition. Models with integer values of sigma = 1 or 2 do not have any transitions. Models with sigma = radical3 and sigma = 3 have weak first-order fluid-solid transitions while those with sigma = 2 radical2, 2 radical3, and 3 radical2 have strong fluid-solid transitions. Pressure, chemical potential, and density are reported for all models and compared to the results for the continuum, theoretical predictions, and prior simulations when available.

Self-referential Monte Carlo method for calculating the free energy of crystalline solids

A self-referential Monte Carlo method is described for calculating the free energy of crystalline solids. All Monte Carlo methods for the free energy of classical crystalline solids calculate the free-energy difference between a state whose free energy can be calculated relatively easily and the state of interest. Previously published methods employ either a simple model crystal, such as the Einstein crystal, or a fluid as the reference state. The self-referential method employs a radically different reference state; it is the crystalline solid of interest but with a different number of unit cells. So it calculates the free-energy difference between two crystals, differing only in their size. The aim of this work is to demonstrate this approach by application to some simple systems, namely, the face centered cubic hard sphere and Lennard-Jones crystals. However, it can potentially be applied to arbitrary crystals in both bulk and confined environments, and ultimately it could also be very efficient.

The melting temperature of the most common models of water

The melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E,TIP3P,TIP4P, TIP4P/Ew, and TIP5P) is obtained from computer simulations at p = 1 bar. Since the melting temperature of ice I(h) for the TIP4P model is now known [E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phys. Rev. Lett. 92, 255701 (2004)], it is possible to use the Gibbs-Duhem methodology [D. Kofke, J. Chem. Phys. 98, 4149 (1993)] to evaluate the melting temperature of ice I(h) for other potential models of water. We have found that the melting temperatures of ice I(h) for SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P models are T = 190 K, 215 K, 146 K, 232 K, 245 K, and 274 K, respectively. The relative stability of ice I(h) with respect to ice II for these models has also been considered. It turns out that for SPC, SPC/E, TIP3P, and TIP5P the stable phase at the normal melting point is ice II (so that ice I(h) is not a thermodynamically stable phase for these models). For TIP4P and TIP4P/Ew, ice I(h) is the stable solid phase at the standard melting point. The location of the negative charge along the H-O-H bisector appears as a critical factor in the determination of the relative stability between the I(h) and II ice forms. The methodology proposed in this paper can be used to investigate the effect upon a coexistence line due to a change in the potential parameters.

A study of the phase behavior of a simple model of chiral molecules and enantiomeric mixtures

We present a study of the solid-fluid and solid-solid phase equilibrium for molecular models representative of chiral molecules and enantiomeric mixtures. The models consist of four hard sphere interaction sites of different diameters in a tetrahedral arrangement with the fifth hard sphere interaction site at the center of the tetrahedron. The volumetric properties and free energies of the pure enantiomers and binary mixtures were calculated in both fluid and solid phases using isobaric Monte Carlo simulations. The models exhibit essentially ideal solution behavior in the fluid phase with little chiral discrimination. In the solid phase the effects of chirality are much greater. Solid-fluid phase behavior involving the pure enantiomer solids and also racemic compounds was calculated. The calculations indicate that, depending on the relative sizes of the hard sphere interaction sites, packing effects alone can be sufficient to stabilize a racemic compound with respect to the pure enantiomer solids.

Solid–liquid phase coexistence of the Lennard-Jones system through phase-switch Monte Carlo simulation

The phase-switch Monte Carlo method of Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] is extended to enable calculation of solid-liquid phase coexistence for soft potentials. The method directly accesses coexistence information about a system while avoiding simulation of the interfacial region. Order parameters are introduced that allow one to define a path that connects liquid and crystalline phases. Transition matrix methods are employed to bias the sampling such that both phases are sampled in a rapid and efficient manner. Coexistence properties are determined through an analysis of specific volume probability distributions, which are generated naturally during a biased simulation. The approach is demonstrated with the Lennard-Jones system. Finite-size effects are examined and compared to those for the hard sphere system. In addition, two techniques are considered for accounting for long-range interactions. The methodology presented here is general and therefore provides a basis for its application to other soft systems.

Solid–fluid and solid–solid phase equilibrium in a model of n-alkane mixtures

Solid-fluid and solid-solid phase equilibrium for binary mixtures of hard sphere chains modeling n-hexane, n-heptane, and n-octane has been calculated using Monte Carlo computer simulations. Thermodynamic integration was used to calculate the Gibbs free energy and chemical potentials in the solid and fluid phases from pure component reference values. A multiple stage free energy perturbation method was used to calculate the composition derivative of the Gibbs free energy. Equation of state and free energy data for the fluid phase indicate ideal solution behavior. Nonideality is much more significant in the solid phase with only partial solubility of shorter chains in the longer chains and essentially no solubility at the other end of the composition range. The miscibility decreases with increasing chain length difference between the components. For the model of n-hexane/n-octane mixtures solid--solid phase separation has been observed directly in some of the simulations, with the components segregating between the layers of the solid structure. The behavior is similar to that seen in some binary n-alkane mixtures with longer chain lengths but comparable chain length ratios between the components. Such phase separation, although indicated thermodynamically, is not seen directly in the simulations of the n-heptane/n-octane mixture due to the difference in the pure component crystal structures.