Communication: Direct determination of triple-point coexistence through cell model simulation.

In simulations of fluid-solid coexistence, the solid phase is modeled as a constrained system of Wigner-Seitz cells with one particle per cell. This model, commonly referred to as the constrained cell model, is a limiting case of a more general cell model, which is formed by considering a homogeneous external field that controls the number of particles per cell and, hence, the relative stability of the solid against the fluid phase. The generalized cell model provides a link that connects the disordered, fluid phase with the ordered, solid phase. In the present work, the phase diagram of this model is investigated through multicanonical simulations at constant pressure and histogram reweighting techniques for a system of 256 Lennard-Jones particles. The simulation data are used to obtain an estimate of the triple point of the Lennard-Jones system. The triple-point pressure is found to be higher compared to previous work. The likely explanation for this discrepancy is the highly compressible nature of the gas phase.

Simulation of phase boundaries using constrained cell models.

Despite impressive advances, precise simulation of fluid-fluid and fluid-solid phase transitions still remains a challenging task. The present work focuses on the determination of the phase diagram of a system of particles that interact through a pair potential, φ(r), which is of the form φ(r) = 4ε[(σ/r)(2n) - (σ/r)(n)] with n = 12. The vapor-liquid phase diagram of this model is established from constant-pressure simulations and flat-histogram techniques. The properties of the solid phase are obtained from constant-pressure simulations using constrained cell models. In the constrained cell model, the simulation volume is divided into Wigner-Seitz cells and each particle is confined to moving in a single cell. The constrained cell model is a limiting case of a more general cell model which is constructed by adding a homogeneous external field that controls the relative stability of the fluid and the solid phase. Fluid-solid coexistence at a reduced temperature of 2 is established from constant-pressure simulations of the generalized cell model. The previous fluid-solid coexistence point is used as a reference point in the determination of the fluid-solid phase boundary through a thermodynamic integration type of technique based on histogram reweighting. Since the attractive interaction is of short range, the vapor-liquid transition is metastable against crystallization. In the present work, the phase diagram of the corresponding constrained cell model is also determined. The latter is found to contain a stable vapor-liquid critical point and a triple point.

Communication: phase transitions, criticality, and three-phase coexistence in constrained cell models.

In simulation studies of fluid-solid transitions, the solid phase is usually modeled as a constrained system in which each particle is confined to move in a single Wigner-Seitz cell. The constrained cell model has been used in the determination of fluid-solid coexistence via thermodynamic integration and other techniques. In the present work, the phase diagram of such a constrained system of Lennard-Jones particles is determined from constant-pressure simulations. The pressure-density isotherms exhibit inflection points which are interpreted as the mechanical stability limit of the solid phase. The phase diagram of the constrained system contains a critical and a triple point. The temperature and pressure at the critical and the triple point are both higher than those of the unconstrained system due to the reduction in the entropy caused by the single occupancy constraint.

Simulation of fluid-solid coexistence via thermodynamic integration using a modified cell model.

Despite recent advances, precise simulation of fluid-solid transitions still remains a challenging task. Thermodynamic integration techniques are the simplest methods to study fluid-solid coexistence. These methods are based on the calculation of the free energies of the fluid and the solid phases, starting from a state of known free energy which is usually an ideal-gas state. Despite their simplicity, the main drawback of thermodynamic integration techniques is the large number of states that must be simulated. In the present work, a thermodynamic integration technique, which reduces the number of simulated states, is proposed and tested on a system of particles interacting via an inverse twelfth-power potential energy function. The simulations are implemented at constant pressure and the solid phase is modeled according to the constrained cell model of Hoover and Ree. The fluid and the solid phases are linked together by performing constant-pressure simulations of a modified cell model. The modified cell model, which was originally proposed by Hoover and Ree, facilitates transitions between the fluid and the solid phase by tuning a homogeneous external field. This model is simulated on a constant-pressure path for a series of progressively increasing values of the field, thus allowing for direct determination of the free energy difference between the fluid and the solid phase via histogram reweighting. The size-dependent results are analyzed using histogram reweighting and finite-size scaling techniques. The scaling analysis is based on studying the size-dependent behavior of the second- and higher-order derivatives of the free energy as well as the dimensionless moment ratios of the order parameter. The results clearly demonstrate the importance of accounting for size effects in simulation studies of fluid-solid transitions.

Precise simulation of the freezing transition of supercritical Lennard-Jones.

The fluid-solid transition of the Lennard-Jones model is analyzed along a supercritical isotherm. The analysis is implemented via a simulation method which is based on a modification of the constrained cell model of Hoover and Ree. In the context of hard-sphere freezing, Hoover and Ree simulated the solid phase using a constrained cell model in which each particle is confined within its own Wigner-Seitz cell. Hoover and Ree also proposed a modified cell model by considering the effect of an external field of variable strength. High-field values favor configurations with a single particle per Wigner-Seitz cell and thus stabilize the solid phase. In previous work, a simulation method for freezing transitions, based on constant-pressure simulations of the modified cell model, was developed and tested on a system of hard spheres. In the present work, this method is used to determine the freezing transition of a Lennard-Jones model system on a supercritical isotherm at a reduced temperature of 2. As in the case of hard spheres, constant-pressure simulations of the fully occupied constrained cell model of a system of Lennard-Jones particles indicate a point of mechanical instability at a density which is approximately 70% of the density at close packing. Furthermore, constant-pressure simulations of the modified cell model indicate that as the strength of the field is reduced, the transition from the solid to the fluid is continuous below the mechanical instability point and discontinuous above. The fluid-solid transition of the Lennard-Jones system is obtained by analyzing the field-induced fluid-solid transition of the modified cell model in the high-pressure, zero-field limit. The simulations are implemented under constant pressure using tempering and histogram reweighting techniques. The coexistence pressure and densities are determined through finite-size scaling techniques for first-order phase transitions which are based on analyzing the size-dependent behavior of susceptibilities and dimensionless moment ratios of the order parameter.

A Monte Carlo study of the freezing transition of hard spheres.

A simulation method for fluid-solid transitions, which is based on a modification of the constrained cell model of Hoover and Ree, is developed and tested on a system of hard spheres. In the fully occupied constrained cell model, each particle is confined in its own Wigner-Seitz cell. Constant-pressure simulations of the constrained cell model for a system of hard spheres indicate a point of mechanical instability at a density which is about 64% of the density at the close packed limit. Below that point, the solid is mechanically unstable since without the confinement imposed by the cell walls it will disintegrate to a disordered, fluid-like phase. Hoover and Ree proposed a modified cell model by introducing an external field of variable strength. High values of the external field variable favor configurations with one particle per cell and thus stabilize the solid phase. In this work, the modified cell model of a hard-sphere system is simulated under constant-pressure conditions using tempering and histogram reweighting techniques. The simulations indicate that as the strength of the field is reduced, the transition from the solid to the fluid phase is continuous below the mechanical instability point and discontinuous above. The fluid-solid transition of the hard-sphere system is determined by analyzing the field-induced fluid-solid transition of the modified cell model in the limit in which the external field vanishes. The coexistence pressure and densities are obtained through finite-size scaling techniques and are in good accord with previous estimates.

Communication: A simple method for simulation of freezing transitions.

Despite recent advances, precise simulation of freezing transitions continues to be a challenging task. In this work, a simulation method for fluid-solid transitions is developed. The method is based on a modification of the constrained cell model which was proposed by Hoover and Ree [J. Chem. Phys. 47, 4873 (1967)]. In the constrained cell model, each particle is confined in a single Wigner-Seitz cell. Hoover and Ree pointed out that the fluid and solid phases can be linked together by adding an external field of variable strength. High values of the external field favor single occupancy configurations and thus stabilize the solid phase. In the present work, the modified cell model is simulated in the constant-pressure ensemble using tempering and histogram reweighting techniques. Simulation results on a system of hard spheres indicate that as the strength of the external field is reduced, the transition from solid to fluid is continuous at low and intermediate pressures and discontinuous at high pressures. Fluid-solid coexistence for the hard-sphere model is established by analyzing the phase transition of the modified model in the limit in which the external field vanishes. The coexistence pressure and densities are in excellent agreement with current state-of-the-art techniques.

Communication: Tracing phase boundaries via molecular simulation: an alternative to the Gibbs-Duhem integration method.

Precise simulation of phase transitions is crucial for colloid/protein crystallization for which fluid-fluid demixing may be metastable against solidification. In the Gibbs-Duhem integration method, the two coexisting phases are simulated separately, usually at constant-pressure, and the phase boundary is established iteratively via numerical integration of the Clapeyron equation. In this work, it is shown that the phase boundary can also be reproduced in a way that avoids integration of Clapeyron equations. The two phases are simulated independently via tempering techniques and the simulation data are analyzed according to histogram reweighting. The main output of this analysis is the density of states which is used to calculate the free energies of both phases and to determine phase coexistence. This procedure is used to obtain the phase diagram of a square-well model with interaction range 1.15σ, where σ is the particle diameter. The phase boundaries can be estimated with the minimum number of simulations. In particular, very few simulations are required for the solid phase since its properties vary little with temperature.

Spatial updating in the great grand canonical ensemble.

In spatial updating grand canonical Monte Carlo, particle transfers are implemented by examining the local environment around a point in space. In the present work, these algorithms are extended to very high densities by allowing the volume to fluctuate, thus forming a great grand canonical ensemble. Since fluctuations are unbounded, a constraint must be imposed. The constrained ensemble may be viewed as a superposition of either constant-pressure or grand canonical ensembles. Each simulation of the constrained ensemble requires a set of weights that must be determined iteratively. The outcome of a single simulation is the density of states in terms of all its independent variables. Since all extensive variables fluctuate, it is also possible to estimate absolute free energies and entropies from a single simulation. The method is tested on a system of hard spheres and the transition from the fluid to a face-centered cubic crystal is located with high precision.

Parallel canonical Monte Carlo simulations through sequential updating of particles.

In canonical Monte Carlo simulations, sequential updating of particles is equivalent to random updating due to particle indistinguishability. In contrast, in grand canonical Monte Carlo simulations, sequential implementation of the particle transfer steps in a dense grid of distinct points in space improves both the serial and the parallel efficiency of the simulation. The main advantage of sequential updating in parallel canonical Monte Carlo simulations is the reduction in interprocessor communication, which is usually a slow process. In this work, we propose a parallelization method for canonical Monte Carlo simulations via domain decomposition techniques and sequential updating of particles. Each domain is further divided into a middle and two outer sections. Information exchange is required after the completion of the updating of the outer regions. During the updating of the middle section, communication does not occur unless a particle moves out of this section. Results on two- and three-dimensional Lennard-Jones fluids indicate a nearly perfect improvement in parallel efficiency for large systems.

Spatial updating grand canonical Monte Carlo algorithms for fluid simulation: generalization to continuous potentials and parallel implementation.

Spatial updating grand canonical Monte Carlo algorithms are generalizations of random and sequential updating algorithms for lattice systems to continuum fluid models. The elementary steps, insertions or removals, are constructed by generating points in space either at random (random updating) or in a prescribed order (sequential updating). These algorithms have previously been developed only for systems of impenetrable spheres for which no particle overlap occurs. In this work, spatial updating grand canonical algorithms are generalized to continuous, soft-core potentials to account for overlapping configurations. Results on two- and three-dimensional Lennard-Jones fluids indicate that spatial updating grand canonical algorithms, both random and sequential, converge faster than standard grand canonical algorithms. Spatial algorithms based on sequential updating not only exhibit the fastest convergence but also are ideal for parallel implementation due to the absence of strict detailed balance and the nature of the updating that minimizes interprocessor communication. Parallel simulation results for three-dimensional Lennard-Jones fluids show a substantial reduction of simulation time for systems of moderate and large size. The efficiency improvement by parallel processing through domain decomposition is always in addition to the efficiency improvement by sequential updating.

Acceleration of Monte Carlo simulations through spatial updating in the grand canonical ensemble.

A new grand canonical Monte Carlo algorithm for continuum fluid models is proposed. The method is based on a generalization of sequential Monte Carlo algorithms for lattice gas systems. The elementary moves, particle insertions and removals, are constructed by analogy with those of a lattice gas. The updating is implemented by selecting points in space (spatial updating) either at random or in a definitive order (sequential). The type of move, insertion or removal, is deduced based on the local environment of the selected points. Results on two-dimensional square-well fluids indicate that the sequential version of the proposed algorithm converges faster than standard grand canonical algorithms for continuum fluids. Due to the nature of the updating, additional reduction of simulation time may be achieved by parallel implementation through domain decomposition.

Parallel Markov chain Monte Carlo simulations.

With strict detailed balance, parallel Monte Carlo simulation through domain decomposition cannot be validated with conventional Markov chain theory, which describes an intrinsically serial stochastic process. In this work, the parallel version of Markov chain theory and its role in accelerating Monte Carlo simulations via cluster computing is explored. It is shown that sequential updating is the key to improving efficiency in parallel simulations through domain decomposition. A parallel scheme is proposed to reduce interprocessor communication or synchronization, which slows down parallel simulation with increasing number of processors. Parallel simulation results for the two-dimensional lattice gas model show substantial reduction of simulation time for systems of moderate and large size.

Simulation of symmetric tricritical behavior in electrolytes.

Despite extensive experimental, theoretical, and simulation efforts, a unified description of ionic phase transitions and criticality has not yet emerged. In this work, we investigate the phase behavior of the restricted primitive model of electrolyte solutions on the simple cubic lattice using grand canonical Monte Carlo simulations and finite-size scaling techniques. The phase diagram of the system is distinctly different from its continuum-space analog. We find order-disorder transitions for reduced temperatures T* < or approximately 0.51, where the ordered structures resemble those of the NaCl crystal. The order-disorder transition is continuous for 0.15 < or approximately T* < or approximately 0.51 and becomes first order at lower temperatures. The line of first-order transitions is a line of three-phase coexistence between a disordered and two ordered phases. The line of continuous, second-order transitions meets this line of triple points at a tricritical point at T* approximately equal to 0.1475. We locate the line of continuous transitions, and the line of triple points using finite-size scaling techniques. The tricritical temperature is estimated by extrapolation of the size-dependent tricritical temperatures obtained from a sixth-order Landau expansion of the free energy. Our calculated phase diagram is in qualitative agreement with mean-field theories.

Scaling fields and pressure mixing in the Widom-Rowlinson model.

We address the issues of scaling fields and of pressure mixing in the penetrable sphere model. This model has an exact symmetry locus from which analytical results may be derived. Based on exact results, we demonstrate that the scaling fields are analytic functions of temperature and chemical potential only. We conclude that there is no pressure mixing in this model. Our findings are in accord with numerical simulations for the same model.

Acceleration of Markov chain Monte Carlo simulations through sequential updating.

Strict detailed balance is not necessary for Markov chain Monte Carlo simulations to converge to the correct equilibrium distribution. In this work, we propose a new algorithm which only satisfies the weaker balance condition, and it is shown analytically to have better mobility over the phase space than the Metropolis algorithm satisfying strict detailed balance. The new algorithm employs sequential updating and yields better sampling statistics than the Metropolis algorithm with random updating. We illustrate the efficiency of the new algorithm on the two-dimensional Ising model. The algorithm is shown to identify the correct equilibrium distribution and to converge faster than the Metropolis algorithm with strict detailed balance. The main advantages of the new algorithm are its simplicity and the feasibility of parallel implementation through domain decomposition.

Monte Carlo simulations of rigid biopolymer growth processes.

Rigid biopolymers, such as actin filaments, microtubules, and intermediate filaments, are vital components of the cytoskeleton and the cellular environment. Understanding biopolymer growth dynamics is essential for the description of the mechanisms and principles of cellular functions. These biopolymers are composed of N parallel protofilaments which are aligned with arbitrary but fixed relative displacements, thus giving rise to complex end structures. We have investigated rigid biopolymer growth processes by Monte Carlo simulations by taking into account the effects of such "end" properties and lateral interactions. Our simulations reproduce analytical results for the case of N = 2, which is biologically relevant for actin filaments. For the case of N = 13, which applies to microtubules, the simulations produced results qualitatively similar to the N = 2 case. The simulation results indicate that polymerization events are evenly distributed among the N protofilaments, which imply that both end-structure effects and lateral interactions are significant. The effect of different splittings in activation energy has been investigated for the case of N = 2. The effects of activation energy coefficients on the specific polymerization and depolymerization processes were found to be unsubstantial. By expanding the model, we have also obtained a force-velocity relationship of microtubules as observed in experiments. In addition, a range of lateral free-energy parameters was found that yields growth velocities in accordance with experimental observations and previous simulation estimates for the case of N = 13.

Monte carlo study of coulombic criticality in polyelectrolytes.

The role of charges in determining the water solubility of polyelectrolytes, a question of considerable relevance to biology, is currently unresolved. We use computer simulations to study the purely Coulombic phase separation of flexible polyelectrolytes with monovalent counterions in an athermal solvent. In agreement with recent theories we find that the critical temperature for this transition increases with chain length, but that the critical density remains unchanged. We therefore stress that the phase behavior of polyelectrolytes is qualitatively different from uncharged polymers, where the critical density decreases towards zero for long chains.

Asymmetric fluid criticality. I. Scaling with pressure mixing.

The thermodynamic behavior of a fluid near a vapor-liquid and, hence, asymmetric critical point is discussed within a general "complete" scaling theory incorporating pressure mixing in the nonlinear scaling fields as well as corrections to scaling. This theory allows for a Yang-Yang anomaly in which mu(")(sigma)(T), the second temperature derivative of the chemical potential along the phase boundary, diverges like the specific heat when T-->T(c); it also generates a leading singular term, /t/(2beta), in the coexistence curve diameter, where t[triple bond](T-T(c))/T(c). The behavior of various special loci, such as the critical isochore, the critical isotherm, the k-inflection loci, on which chi((k))[triple bond]chi(rho,T)/rho(k) (with chi=rho(2)k(B)TK(T)) and C((k))(V)[triple bond]C(V)(rho,T)/rho(k) are maximal at fixed T, is carefully elucidated. These results are useful for analyzing simulations and experiments, since particular, nonuniversal values of k specify loci that approach the critical density most rapidly and reflect the pressure-mixing coefficient. Concrete illustrations are presented for the hard-core square-well fluid and for the restricted primitive model electrolyte. For comparison, a discussion of the classical (or Landau) theory is presented briefly and various interesting loci are determined explicitly and illustrated quantitatively for a van der Waals fluid.

Precise simulation of criticality in asymmetric fluids.

Extensive grand canonical Monte Carlo simulations have been performed for the hard-core square-well fluid with interaction range b=1.5 sigma. The critical exponent for the correlation length has been estimated in an unbiased fashion as nu=0.63+/-0.03 via finite-size extrapolations of the extrema of properties measured along specially constructed, asymptotically critical loci that represent pseudosymmetry axes. The subsequent location of the critical point achieves a precision of five parts in 10(4) for Tc and about 0.3% for the critical density rhoc. The effective exponents gamma+(eff) and beta(eff) indicate Ising-type critical-point values to within 2% and 5.6%, respectively, convincingly distinguishing the universality class from the "nearby" XY and n=0 (self-avoiding walk) classes. Simulations of the heat capacity CV(T,rho) and d2psigma/dT2, where psigma is the vapor pressure below Tc, suggest a negative but small Yang-Yang anomaly, i.e., a specific-heat-like divergence in the corresponding chemical potential derivative (d2 musigma/dT2) that requires a revision of the standard asymptotic scaling description of asymmetric fluids.