Canonical Ensemble Simulations Research Articles

Self-Adapting Fixed-End-Point Configurational-Bias Monte Carlo Method for the Regrowth of Interior Segments of Chain Molecules with Strong Intramolecular Interactions

An extension to the configurational-bias Monte Carlo method is presented which allows for the efficient conformational sampling of the interior segments of chain molecules whose interactions include strong bonded terms (governing bond stretching, bond angle bending, and dihedral angle rotation). The ability to regrow interior segments overcomes the limitations of conventional configurational-bias methods (where the regrowth is always directed to a free chain end) and now allows for the simulation of chain systems with low concentrations of chain ends, that is, higher molecular weights, networks, or cyclic structures. As previously proposed by Dijkstra et al. [J. Chem. Phys. 1994, 101, 3179] for lattice polymers and by Vendruscolo [J. Chem. Phys. 1997, 106, 2970] for freely jointed polymers, an additional biasing (closing) probability is used that guides the bead-by-bead configurational-bias regrowth of interior segments toward its desired fixed target. However, while the previous methods are limited to chain models for which the number of random walks that lead to closure is known or which rely on simpler and less efficient geometric considerations, the algorithm presented here allows for the simulation of chain molecules using force fields of arbitrary complexity for which the closing probability is not known a priori. It is important to note that the additional biasing probability used to guide the move does not necessarily have to be the true closing probability but that a good approximation thereof is essential to improve the sampling efficiency. To this extent, we obtain an intial guess of the biasing probability from a short presimulation or an earlier simulation of a related system or simply use a uniform biasing probability. A self-adapting scheme is then used to optimize the biasing probability during the course of the simulation for the system of interest. In addition to the conformational sampling of interior segments, the new algorithm also enables efficient particle insertions and removals of cyclic molecules (of moderate length) and thereby opens the door to simulations in the grand canonical and Gibbs ensembles. Simulation results are presented for linear, branched, and cyclic alkanes using the transferable potentials for phase equilibria (TraPPE) force field.

Renormalized Hamiltonian for a peptide chain: Digitalizing the protein folding problem

A renormalized Hamiltonian for a flexible peptide chain is derived to generate the long-time limit dynamics compatible with a coarsening of torsional conformation space. The renormalization procedure is tailored taking into account the coarse graining imposed by the backbone torsional constraints due to the local steric hindrance and the local backbone-side-group interactions. Thus, the torsional degrees of freedom for each residue are resolved modulo basins of attraction in its so-called Ramachandran map. This Ramachandran renormalization (RR) procedure is implemented so that the chain is energetically driven to form contact patterns as their respective collective topological constraints are fulfilled within the coarse description. In this way, the torsional dynamics are digitalized and become codified as an evolving pattern in a binary matrix. Each accepted Monte Carlo step in a canonical ensemble simulation is correlated with the real mean first passage time it takes to reach the destination coarse topological state. This real-time correlation enables us to test the RR dynamics by comparison with experimentally probed kinetic bottlenecks along the dominant folding pathway. Such intermediates are scarcely populated at any given time, but they determine the kinetic funnel leading to the active structure. This landscape region is reached through kinetically controlled steps needed to overcome the conformational entropy of the random coil. The results are specialized for the bovine pancreatic trypsin inhibitor, corroborating the validity of our method.

Interface and Surface Properties of Short Polymers in Solution: Monte Carlo Simulations and Self-Consistent Field Theory

We investigate the structure and thermodynamics of inhomogeneous polymer solutions in the framework of a coarse-grained off-lattice model. Properties of the liquidvapor interface and the packing of the solution in contact with an attractive wall are considered. The wall interacts with the monomers via a long-range interaction. We employ Monte Carlo simulations in the grand canonical ensemble and self-consistent field calculations to study segment profiles, excess free energies, and the wetting behavior. The self-consistent field calculations incorporate the detailed chain structure via a partial enumeration scheme and use a weighteddensity functional ansatz for the monomeric interactions. Quantitative agreement between the Monte Carlo simulations and the self-consistent field calculations is achieved for the conformational arrangements at the wall. Only qualitative agreement is obtained for the excess free energies of interfaces and surfaces. Possible sources for the quantitative deviations are discussed. Using the Young equation, we accurately locate the first-order wetting transition. In the Monte Carlo simulations and the self-consistent field calculations, we find a strong first-order wetting transition and locate the concomitant prewetting line. The behavior of the prewetting line close to the bulk coexistence yields a compatible estimate for locating the wetting transition. The location of the prewetting critical point is estimated in the simulations.

Adiabatic Nuclear and Electronic Sampling Monte Carlo Simulations in the Gibbs Ensemble: Application to Polarizable Force Fields for Water

The adiabatic nuclear and electronic sampling Monte Carlo algorithm (ANES-MC) is extended to simulations in the Gibbs ensemble. Whereas the maximum displacements used for translational, rotational, and volume trial moves can be adjusted to foster efficient sampling in the adiabatic limit, the transfer (swap) of particles always causes a major disturbance of the electronic structures of the two phases (supplying and receiving the particle). To reequilibrate the electronic structures requires additional sampling of the electronic degrees of freedom. A simple, distance-dependent criterion for the preferential selection of the electronic degrees of freedom, for which a move is to be attempted, is shown to improve the efficiency of the particle swap move. The ANES-MC algorithm is applied to the polarizable simple point charge-fluctuating charge (SPC−FQ) and transferable intermolecular potential 4 point-fluctuating charge (TIP4P−FQ) models proposed by Rick et al. (J. Chem. Phys. 1994, 101, 6141). For both models simulations were performed using the standard constraint on the neutrality of individual molecules. In addition, for the SPC−FQ model the use of a constraint on the neutrality of an entire phase was investigated, which allows for intermolecular charge transfer. Simulations in the Gibbs ensemble were performed to calculate the vapor−liquid coexistence curves from 323 to 523 K, whereas simulations in the grand canonical ensemble were performed for the near-critical region. Dielectric constants at different state points were calculated from canonical ensemble simulations. Neither the SPC−FQ nor the TIP4P−FQ force fields give a satisfactory description of the vapor−liquid equilibria. In particular, the critical temperature is greatly underestimated by both models. Although intermolecular charge transfer has only a very small influence on the internal energy and the radial distribution functions at ambient conditions and along the coexistence curve, it increases the dielectric constant by approximately 30%.

Comment on “Molecular dynamics simulations in the grand canonical ensemble: Formulation of a bias potential for umbrella sampling” [J. Chem. Phys. 110, 8295 (1999)

The cavity-biased grand canonical ensemble Monte Carlo simulation method has been compared with a particular implementation of the grand canonical ensemble molecular dynamics method by examining the fluctuations during simulations of equivalent length. The Monte Carlo run exhibited an order of magnitude faster fluctuation than the molecular dynamics run.

Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling

A grand canonical ensemble molecular dynamics (GMD) simulation method has been adapted to examine the thermodynamics of clay-mineral hydration. In the GMD method, the number of water molecules in the system is treated as a continuous variable for which an equation of motion is established. Fluctuations in the water content at constant chemical potential are investigated using trajectories of this particle number variable. A bias potential may be used to modify the free energy contour along the particle number coordinate. This catalyzes particle fluctuations and greatly improves simulation convergence. Adaptation of the GMD method to treat hydrated clay minerals included the introduction of a local-control technique that fixes the water chemical potential in the clay interlayer region. In addition, a bias-potential feedback algorithm was implemented to improve particle fluctuation efficiency. Information pertaining to the free energy contour, generated during the course of the simulation, was used periodically to enhance the bias potential. This allowed for the utilization of a single input bias potential under a broad range of simulation conditions. The method was used to investigate swelling of a cesium–montmorillonite clay. Measured disjoining pressures showed oscillations that are indicative of crystalline-swelling phase transitions. Integration of the disjoining pressures yielded a swelling free energy profile with distinct free-energy minima for the one- and two-layer hydrates. The results may be compared qualitatively with both clay swelling and surface force apparatus experiments, and with previous simulation studies of simple fluids in slit pores.

Molecular dynamics simulations in the grand canonical ensemble: Formulation of a bias potential for umbrella sampling

An extended Hamiltonian technique for performing grand canonical ensemble molecular dynamics simulations has been reformulated to include umbrella sampling, thus improving the efficiency of particle creation and annihilation processes. This was accomplished through incorporation of a bias potential in the Hamiltonian that modifies the free energy contour between integer particle number states. The extended Hamiltonian includes a continuous particle number variable that is the sum of the integer particle number and a coupling parameter, the latter being used in a scaling function for a fractional molecule’s interactions with the rest of the system. Equations of motion for the coupling parameter, derived from the Hamiltonian, integrate to yield density fluctuations at constant chemical potential. This new method may be adapted to a wide range of potential scaling functions. The technique was applied to calculations of extended simple point charge water density versus chemical potential, using both linear and nonlinear scaling. For each scaling function, a bias potential was constructed using a thermodynamic integration technique. Grand canonical ensemble simulations then yielded results in agreement with independent calculations.

On the Scaling of the Critical Solution Temperature of Binary Polymer Blends with Chain Length

Monte Carlo simulations have been performed to examine the scaling with chain length of the upper critical solution temperature (Tc) of binary, symmetric blends of polymers. Critical parameters were obtained by histogram reweighting analysis of semigrand canonical ensemble simulations near the critical temperature. For several continuum space polymeric models, it is found that Tc/ε scales approximately as either ρc(A + B N) or ρcN/(B-1 + A/ ), where N is chain length, ρc is the density at Tc, ε is a chemical mismatch parameter, and A and B are constants that depend on the specific characteristics of the model. Constant A provides a measure of the correction to the Flory−Huggins' prediction Tc ∼ N. The effect of such correction term becomes unimportant for large N, ε, and temperature. Excellent agreement is found, however, with the mean-field prediction that χE ∼ 2/N (for all systems studied), where χE is a continuum space, enthalpic “chi” parameter which takes into account the variations of local structure of the fluid for different chain lengths at the critical point.

Phase behavior of the restricted primitive model and square-well fluids from Monte Carlo simulations in the grand canonical ensemble

Coexistence curves of square-well fluids with variable interaction width and of the restricted primitive model for ionic solutions have been investigated by means of grand canonical Monte Carlo simulations aided by histogram reweighting and multicanonical sampling techniques. It is demonstrated that this approach results in efficient data collection. The shape of the coexistence curve of the square-well fluid with short potential range is nearly cubic. In contrast, for a system with a longer potential range, the coexistence curve closely resembles a parabola, except near the critical point. The critical compressibility factor for the square-well fluids increases with increasing range. The critical behavior of the restricted primitive model was found to be consistent with the Ising universality class. The critical temperature was obtained as Tc=0.0490±0.0003 and the critical density ρc=0.070±0.005, both in reduced units. The critical temperature estimate is consistent with the recent calculation of Caillol et al. [J. Chem. Phys. 107, 1565 (1997)] on a hypersphere, while the critical density is slightly lower. Other previous simulations have overestimated the critical temperature of this ionic fluid due to their failure to account for finite-size effects in the critical region. The critical compressibility factor (Zc=Pc/ρcTc) for the ionic fluid was obtained as Zc=0.024±0.004, an order of magnitude lower than for nonionic fluids.

Critical point and phase behavior of the pure fluid and a Lennard-Jones mixture

Monte Carlo simulations in the grand canonical ensemble were used to obtain liquid-vapor coexistence curves and critical points of the pure fluid and a binary mixture of Lennard-Jones particles. Critical parameters were obtained from mixed-field finite-size scaling analysis and subcritical coexistence data from histogram reweighting methods. The critical parameters of the untruncated Lennard-Jones potential were obtained as Tc*=1.3120±0.0007, ρc*=0.316±0.001 and pc*=0.1279±0.0006. Our results for the critical temperature and pressure are not in agreement with the recent study of Caillol [J. Chem. Phys. 109, 4885 (1998)] on a four-dimensional hypersphere. Mixture parameters were ε1=2ε2 and σ1=σ2, with Lorentz–Berthelot combining rules for the unlike-pair interactions. We determined the critical point at T*=1.0 and pressure-composition diagrams at three temperatures. Our results have much smaller statistical uncertainties relative to comparable Gibbs ensemble simulations.

Hard sphere mixtures near a hard wall

Hard spheres with ratios of diameters of 3:5 and 1:3 and at various densities are studied by means of grand canonical ensemble (GCE) simulations, first- and second-order Percus—Yevick theory, scaled particle theory (SPT) and density functional theory (DFT). Generally, but not always, the density profiles of the first-order Percus—Yevick (PY1) theory results are smaller at contact than the GCE simulation results. The second-order Percus—Yevick (PY2) and DFT results are in much better agreement with the simulations. PY2 theory seems to be more accurate than DFT but is much more demanding computationally.

The chemical potential of Lennard-Jones associating fluids from osmotic Monte Carlo simulations

The chemical potential for a two-component Lennard-Jones fluid with associative interaction between opposite species promoting the formation of dimers is calculated using osmotic Monte Carlo (OMC) canonical ensemble simulations. Grand canonical Monte Carlo simulations also are performed to verify the accuracy of the OMC approach. The data from both methods agree very well for thermodynamic states with different degrees of dimerization. It follows that the OMC is a promising approach for the determination of the thermodynamics of and equilibria between associating and non-associating fluids and associating fluid mixtures.

Confined thin films of linear and branched alkanes

We report computer simulations in the Grand canonical ensemble of a system of several alkanes between two solid surfaces. We computed the solvation force exerted by the fluid on the plates. The solvation force for linear decane oscillates with distance with a periodicity close to the width of the molecules. The branched alkanes (2-methylundecane and 2-methylheptane) show a similar oscillatory behavior, however the oscillations are decreased and are shifted to the attractive regime. In addition, we computed the liquid-vapour equilibria by using Gibbs ensemble Monte-Carlo simulations of n-pentane confined in a slit of 9, 13, 17 Å. The critical temperature of the liquid-vapour coexistence shifts to lower temperatures upon confining. At a plate separation of 5 Å, no liquid-vapour equilibrium is found.

Non-neutrality in a charged slit pore

Grand canonical Monte Carlo simulations are reported for restricted primitive model ions confined between two parallel charged walls. Single-ion insertion/deletion in the grand canonical ensemble simulation allows individual ions inside the slit pore to equilibrate with an external solution. Due to the confinement effect, the net charge of the ions inside is less than that required to balance the wall charges. The residue charge increases with a narrower slit or a higher charge density on the walls. The ion distribution between the walls is reported and compared with the modified Gouy Chapman results.

Adsorption from binary solution in slit shaped pores

AbstractA theoretical model for adsorption of a binary solution in a slit shaped pore is proposed. The classic lattice model of the solution is used. Adsorption systems exhibiting complete mixing are considered. The theoretical results are compared with Monte Carlo simulations in the grand canonical ensemble. An influence of pore size and non‐ideality of the liquid mixture on adsorption and selectivity is analysed.

Monte Carlo simulation of athermal mesogenic chains: Pure systems, mixtures, and constrained environments

Monte Carlo simulations of mesogenic polymers have been performed for single and binary systems in the bulk and for pure components in constrained environments. An expanded ensemble formalism was employed to facilitate fluctuations of the number of particles in grand canonical, Gibbs, and pseudo-Gibbs ensemble simulations. Phase diagrams were generated for pure systems of athermal 8-mers and 16-mers with varying degrees of flexibility, for a binary mixture of rigid and semiflexible 16-mers (at constant pressure and variable flexibility), and for a bidisperse mixture of rigid chains. Simulations of the partitioning of rigid 12-mer chains between a bulk and a slit pore revealed that, for a certain range of slit gaps, the preferential coplanar orientation of molecules at the proximity of the walls induces the formation of a nematic phase within the pore. A similar induction of a nematic phase was observed for anisotropic molecules dispersed in a stretched, swollen gel. For the bulk systems, our simulation results were compared to predictions of various theories.

Monte Carlo Simulations of Fluid Phase Transitions in Micropores: Hysteresis and Van Der Waals-Type Loops

ABSTRACTWe report the results of the Monte Carlo simulation study of the fluid phase transitions in micropores. The Gibbs, Canonical, and Grand Canonical ensembles were employed to obtain the adsorption isotherms of the Lennard-Jones nitrogen in carbon slit micropores at 77.4 K. It is shown for the first time that depending on the choice of the parameters of the Monte Carlo simulations, such as the number of particles in the simulation cells, the sampling from the Gibbs ensemble may supply either van der Waals-type or hysteresis loops in the phase transition zones. The existence of these loops is confirmed by independent sampling from the Canonical ensemble using Widom method of test particles. The Grand Canonical ensemble simulations always produce hysteresis loops within the phase transition zone. Outside the phase transition zones, sampling from all three ensembles yield essentially identical isotherms.

Monte Carlo simulations of a polyoxyethylene C12E2 lamellar bilayer in water

The polyoxyethylene C12E2/ water lamellar phase has been modelled as an aqueous region sandwiched by two bilayers. The bilayers were represented as ethoxy chains attached to a hydrocarbon continuum. Initially, the bilayer separation distance was fixed at 24 A, which is a typical experimental value for the surfactant weight fraction of 71%. Both canonical and Gibbs ensemble simulations have been performed on the system at different temperatures ranging from 298⋅15 K to 353⋅15 K. In the Gibbs simulations the bilayer separation was allowed to readjust to an equilibrium value. The results are compared and contrasted with previous studies in which single chains (Kong, Y. C., Nicholson, D., Parsonage, N. G., and Thompson, L., 1994, J. chem Soc. Faraday Trans, 90, 2375) and charged OH groups (Cracknell, R. F., Nicholson, D., and Parsonage, N. G., 1992, Molec. Phys., 75, 1023), attached to a hydrocarbon substrate, were simulated. As in the previous work with single chains, there is evidence of water bridging whic...

Grand canonical ensemble Monte Carlo simulation of the dCpG/proflavine crystal hydrate

The grand canonical ensemble Monte Carlo molecular simulation method is used to investigate hydration patterns in the crystal hydrate structure of the dCpG/proflavine intercalated complex. The objective of this study is to show by example that the recently advocated grand canonical ensemble simulation is a computationally efficient method for determining the positions of the hydrating water molecules in protein and nucleic acid structures. A detailed molecular simulation convergence analysis and an analogous comparison of the theoretical results with experiments clearly show that the grand ensemble simulations can be far more advantageous than the comparable canonical ensemble simulations.

Simulations of a confined polymer fluid at constant chemical potential

Monte Carlo simulations are reported of a polymer fluid confined between planar walls using the isotension ensemble, which allows one to perform simulations at constant chemical potential without the necessity of particle exchanges with a bulk polymer fluid. The pressure acting upon the walls and the free energy per unit area of the confined fluid as a function of separation between the walls are obtained, all at constant chemical potential. This is made possible by the use of the so-called free energy difference method. This method appears to be less susceptible to the difficulties that confront the simulation of confined fluids at high density in grand canonical ensemble simulations.