Information theory approach to the isothermal–isobaric ensemble. 1. Monte Carlo simulations with the hard-sphere and square potentials

Monte Carlo (MC) simulations were performed on the isothermal–isobaric partition functions for both argon and methane gas. A newly implemented form was applied to the calculation of the volume for a variety of pressures, from which many potential applications can be derived.     

Monte Carlo simulations of squalane in the Gibbs ensemble

We investigate the liquid–vapour coexistence curve of 2,6,10,15,19,23-hexamethyltetracosane (squalane) near the critical point with a new Lennard–Jones parameter set and compare our results to existing simulation data as well as to recent experimental vapour pressure data. Comparison of the liquid–vapour coexistence curve to previous simulation data reveals that this new force field, which includes tail corrections to the truncation of the non-bonded interactions increases the liquid density. We determine the critical temperature to 829 K and 825 K (with roughly 1% error) for two different system sizes, 72 and 108 molecules, and the critical density to 0.211 g/cm³ and 0.228 g/cm³, respectively. We extrapolate experimental vapour pressure data by use of Antoine's law to the temperature range covered by simulation and yield good agreement between simulation and experiment. We note that the vapour pressure in simulation is essentially governed by the ideal vapour pressure.     

Optimized ensemble Monte Carlo simulations of dense Lennard-Jones fluids

We apply the recently developed adaptive ensemble optimization technique to simulate dense Lennard-Jones fluids and a particle-solvent model by broad-histogram Monte Carlo techniques. Equilibration of the simulated fluid is improved by sampling an optimized histogram in radial coordinates that shifts statistical weight towards the entropic barriers between the shells of the liquid. Interstitial states in the vicinity of these barriers are identified with unprecedented accuracy by sharp signatures in the quickly converging histogram and measurements of the local diffusivity. The radial distribution function and potential of mean force are calculated to high precision.     

Convergence of Monte Carlo simulations of liquid water in the NPT ensemble

Long Monte Carlo simulations of liquid water at 25° and 1 atm have been carried out to study the convergence characteristics of the calculations. The recently reported TIPS2 potential was employed with system sizes of 125 and 216 monomers. The enthalpy, volume and radial distribution functions converge rapidly and show little size dependence. However, the rates of convergence are much slower for the fluctuation properties and are in the order: heat capacity (Cp) ? isothermal compressibility (?) ? coefficient of thermal expansion (α). In fact,the weak coupling of the enthalpy and volume allows only crude estimates for α. In addition, the estimation of statistical error bounds for the thermodynamic properties is analyzed. It is found that very long (2500-3500 k) simulations are needed to yield valid estimates of the errors for the enthalpy, volume and C_p, while the erros for ? and α would still be elusive.     

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.     

Diffusion Monte Carlo simulations on uracil-water using an anisotropic atom-atom potential model

We have developed an anisotropic atom-atom intermolecular potential model for the interaction of uracil with water. The potential consists of a distributed multipole analysis (DMA) model for the electrostatic energy, and a 6-exp potential to represent the repulsion-dispersion term. The repulsion-dispersion potential parameters are adjusted to yield good agreement with accurate ab initio data on the minima and transition states of the uracil-water complex. We have used this potential in diffusion Monte Carlo simulations of uracil-water, uracil-(water)2 and uracil-(water)3. The uracil-water simulations show that the theoretically based potential gives a qualitatively different picture of uracil hydration than that provided by a standard isotropic atom-atom point charge model, which is shown to underestimate the delocalized motion of the water hydrogen atoms. Plots of the vibrational probability density of the hydrogen atoms show the delocalized motion of the water hydrogen atoms that are not involved in hydrogen bonding.     

Calculating Gibbs free energies of transfer from Gibbs ensemble Monte Carlo simulations

The partitioning of the ternary systems n-pentane/n-heptane/(helium or argon) at ambient conditions is investigated using configurational-bias Monte Carlo simulations in the Gibbs ensemble. The results demonstrate that this approach yields very precise partition constants and free energies of transfer. Simulations are carried out to study the dependence of the n-pentane partitioning with respect to the carrier gas, the system size, and the overall solute concentrations. None of the changes of variables, within the ranges used here, has a significant effect on the alkane partitioning. However, chemical potentials calculated via Widom's ghost particle insertions show a strong number dependence for phases containing relatively few molecules of a given type. This problem originates from the fact that the chemical potential is calculated for a concentration of real particles plus one ghost particle that is systematically larger than the equilibrium concentration. A simple correction term is suggested to account for this problem. Received: 13 May 1998 / Accepted: 18 June 1998 / Published online: 4 September 1998     

Monte Carlo simulations of a liquid crystal copolymer in the solid state

The condensed phase of the alternating copolyester of p-hydroxybenzoic acid (HBA) and 2-hydroxy-6-naphthoic acid (HNA) is investigated by studying the room temperature packing arrangement of the copolymer chains. A molecular modeling methodology is employed with a Monte Carlo sampling of the configurational phase space. Realistic poly(HBA-alt-HNA) polymer chains are represented by an explicit atom representation of the HBA/HNA dimers. States are sampled from the NVT ensemble using a sampling scheme consisting of (1) valence and torsional variations, (2) rigid body rotations of the chain about the chain axis, and (3) rigid body translations of the chain. The effect of chain packing on the conformation of chains, as well as the relative intra- and intermolecular orientations of aromatic rings, is investigated. Correlation of chain positioning along the chain axis is dominated by aromatic rings maintaining a center-to-center plane of registry. These layers of aromatic units pack with a preference for edge-to-face orientations in a herringbone-type pattern and have an intermolecular ring angle between the pairs of aromatic rings in the unit cell that is ca. 68°. The aromatic rings, on average, are rotated 38° out from the b–c plane. The phenylene rings of these copolyesters are less restricted in their relative orientation in comparison to the naphthalene rings. Intramolecular orientational probability density distributions indicate a preference for staggering the successive aromatic rings along the chain, with a staggering angle of ca. 66°. © 1998 John Wiley & Sons, Inc. J Polym Sci B: Polym Phys 36: 727–741, 1998     

Monte Carlo simulations of the hydration of substituted benzenes with OPLS potential functions

Intermolecular potential functions have been developed for use in computer simulations of substituted benzenes. Previously reported optimized potentials for liquid simulations (OPLS) for benzene and organic functional groups were merged and tested by computing free energies of hydration for toluene, p-xylene, phenol, anisole, benzonitrile, p-cresol, hydroquinone, and p-dicyanobenzene. The calculations featured Monte Carlo simulations at 25°C and 1 atm with statistical perturbation theory. The average difference between the computed results and experimental data for the absolute free energies of hydration is 0.5 kcal/mol. The AM1-SM2 method is also found to perform well in predicting the free energies of hydration for the substituted benzenes. In addition, the Monte Carlo simulations provided details on the hydration of the substituted benzenes, in particular for the solute–water hydrogen bonding. © 1993 John Wiley & Sons, Inc.     

Monte Carlo simulations of pure liquid substituted benzenes with OPLS potential functions

Intermolecular potential functions have been developed for use in computer simulations of substituted benzenes. Previously reported optimized potentials for liquid simulations (OPLS) for benzene and organic functional groups were merged and tested in Monte Carlo statistical mechanics simulations for the pure liquids of toluene, m-cresol, anisole, aniline, and benzonitrile at 25°C at 1 atm. The merged potential functions yielded acceptable thermodynamic results for the liquids except in the case of aniline, for which the error in the heat of vaporization was 12%. This was remedied by enhancing the polarity of the model to be more consistent with the observed dipole moment of aniline. Overall, the average errors in computed heats of vaporization and densities were then 2 and 1%, respectively. The structures of the liquids were characterized through energy and radial distribution functions. For m-cresol and aniline, the molecules participate in averages of 1.6 and 1.4 hydrogen bonds, respectively. Condensed phase effects on the torsional energies for anisole, m-cresol, and aniline were found to be small; m-cresol has a slightly enhanced tendency to be nonplanar in the liquid than in the gas phase, while anisole shows the opposite pattern. © 1993 John Wiley & Sons, Inc.     

Monte Carlo simulations of fluids whose particles interact with a logarithmic potential

Monte Carlo simulations of a model fluid in which the particles interact via a continuous potential that has a logarithmic divergence at a pair separation of sigma, which we introduced in J. G. Powles et al., Proc. R. Soc. London, Ser. A 455, 3725 (1999), have been carried out. The potential has the form, phi(r)= -epsilon ln(fr), where epsilon sets the energy scale and fr=1-(sigma/r)m. The value of m chosen was 12 but the qualitative trends depend only weakly on the value of m, providing it is greater than 3. The potential is entirely repulsive and has a logarithmic divergence as approximately -ln(r/sigma-1) in the r-->sigma limit. Predictions of the previous paper that the internal energy can be computed at all temperatures using the standard statistical mechanics formula for continuous potentials are verified here. The pressure can be calculated using the usual virial expression for continuous potentials, although there are practical limitations in resolving the increasingly important contribution from the r-->sigma limit at reduced temperatures greater than approximately 5. The mean square force F2 and infinite frequency shear Ginfinity and bulk Kinfinity moduli are only finite for T*=kBT/epsilon<1. The logarithmic fluid's physical properties become increasingly more like that of the hard sphere fluid with increasing temperature, showing a sharp transition in the behavior of the mean square force and infinite frequency elastic constants at T*=1. The logarithmic fluid is shown to exhibit a solid-fluid phase transition.     

Donnan potential of dilute colloidal dispersions: Monte Carlo simulations

Donnan equilibrium between a salt-free colloidal dispersion and an electrolyte solution has been investigated by Monte Carlo simulations. The Donnan potential is directly calculated by considering two compartments separated by a semipermeable membrane. In order to understand the role played by colloid–ion interactions, the influences of colloidal characteristics, including particle size R, intrinsic particle charge Z, counterion valency z_c, and concentration c_p, on Donnan potential Ψ_D and effective charge Z_eff are examined. Our simulations show that the electroneutrality condition is not followed in both compartments and the Donnan potential arises due to the net charge density. The Donnan potential grows by increasing c_p and Z_eff and by decreasing dielectric constant ε_r, i.e., Ψ_DZ_effc_p/ε_r approximately. Note that the effective charge varies with R,Z,c_p,ε_r and z_c as well. When the salt concentration is increased, the net charge density is lowered and thus the Donnan potential decays accordingly. The validity of the classical theory based on the Nernst equation and the electroneutrality assumption is also examined. In general, the simulation results at high-salt condition can be well represented by such mean-field theory.     

Water in nanopores. I. Coexistence curves from Gibbs ensemble Monte Carlo simulations

Coexistence curves of water in cylindrical and slitlike nanopores of different size and water-substrate interaction strength were simulated in the Gibbs ensemble. The two-phase coexistence regions cover a wide range of pore filling level and temperature, including ambient temperature. Five different kinds of two-phase coexistence are observed. A single liquid-vapor coexistence is observed in hydrophobic and moderately hydrophilic pores. Surface transitions split from the main liquid-vapor coexistence region, when the water-substrate interaction becomes comparable or stronger than the water-water pair interaction. In this case prewetting, one and two layering transitions were observed. The critical temperature of the first layering transition decreases with strengthening water-substrate interaction towards the critical temperature expected for two-dimensional systems and is not sensitive to the variation of pore size and shape. Liquid-vapor phase transition in a pore with a wall which is already covered with two water layers is most typical for hydrophilic pores. The critical temperature of this transition is very sensitive to the pore size, in contrast to the liquid-vapor critical temperature in hydrophobic pores. The observed rich phase behavior of water in pores evidences that the knowledge of coexistence curves is of crucial importance for the analysis of experimental results and a prerequiste of meaningful simulations.     

Monte Carlo simulation of liquid hydroxylamine using ab initio intermolecular potential functions

A potential model for intermolecular interactions between hydroxylamine (NH₂OH) molecules based on ab initio quantum mechanical calculations is reviewed and critically assessed by analyzing results from a Monte Carlo simulation of liquid hydroxylamine. The liquid structure is studied in detail using radial, energy, and angular distribution functions, coordination numbers, and their distribution. Results indicate a large first solvation shell (5.3 Å), which contains 13 molecules, out of which only 4 are truly bonded by nonlinear, low-energy hydrogen bonds. These are of either the OH…O or the OH…N type, as NH…O and NH…N linear bonds are considerably suppressed, and no cyclic dimers are found. The dependence of the structural and physical properties on the simulation characteristics has also been investigated.     

Stochastic potential switching algorithm for Monte Carlo simulations of complex systems

This paper describes a new Monte Carlo method based on a novel stochastic potential switching algorithm. This algorithm enables the equilibrium properties of a system with potential V to be computed using a Monte Carlo simulation for a system with a possibly less complex stochastically altered potential V. By proper choices of the stochastic switching and transition probabilities, it is shown that detailed balance can be strictly maintained with respect to the original potential V. The validity of the method is illustrated with a simple one-dimensional example. The method is then generalized to multidimensional systems with any additive potential, providing a framework for the design of more efficient algorithms to simulate complex systems. A near-critical Lennard-Jones fluid with more than 20,000 particles is used to illustrate the method. The new algorithm produced a much smaller dynamic scaling exponent compared to the Metropolis method and improved sampling efficiency by over an order of magnitude.     

Monte Carlo simulations of the exchange kinetics of polymers adsorbed on a solid surface

The exchange kinetics of polymers adsorbing on a solid surface is extensively studied by dynamic Monte Carlo simulations. A model employed simulates a semidilute polymer solution placed in contact with a solid surface that attracts polymer segments by the adsorption interaction (χ_s). The exchange process of polymer chains, between the solution and the adsorbed polymer layer, is examined under various conditions. The exchange kinetics shows two characteristic regimes with increasing chain length. One is the diffusion‐controlled regime found with a small χ_s , and the other the detachment‐controlled regime with a large χ_s . These two regimes are well described by a kinetic theory. Various dynamic quantities show that the diffusion‐controlled regime is not due to sluggish dynamics near the surface, but rather to bulk diffusion of chains. The diffusion‐controlled regime found in this study is considered to appear at the high temperature limit.     

A theoretical study of the hydration of Rb+ by Monte Carlo simulations with refined ab initio-based model potentials

An infinitely diluted aqueous solution of Rb⁺ was studied using ab initio-based model potentials in classical Monte Carlo simulations to describe its structural and thermodynamic features. An existing flexible and polarizable model [Saint-Martin et al. in J Chem Phys 113(24) 10899, 2000] was used for water–water interactions, and the parameters of the Rb⁺–water potential were fitted to reproduce the polarizability of the cation and a sample of ab initio pair interaction energies. It was necessary to calibrate the basis set to be employed as a reference, which resulted in a new determination of the complete basis set (CBS) limit energy of the optimal Rb⁺–OH₂ configuration. Good agreement was found for the values produced by the model with ab initio calculations of three- and four-body nonadditive contributions to the energy, as well as with ab initio and experimental data for the energies, the enthalpies and the geometric parameters of Rb⁺(H₂O) _n clusters, with n = 1,  2,…, 8. Thus validated, the potential was used for simulations of the aqueous solution with three versions of the MCDHO water model; this allowed to assess the relative importance of including flexibility and polarizability in the molecular model. In agreement with experimental data, the Rb⁺–O radial distribution function (RDF) showed three maxima, and hence three hydration shells. The average coordination number was found to be 6.9, with a broad distribution from 4 to 12. The dipole moment of the water molecules in the first hydration shell was tilted to 55° with respect to the ion’s electric field and had a lower value than the average in bulk water; this latter value was recovered at the second shell. The use of the nonpolarizable version of the MCDHO water model resulted in an enhanced alignment to the ion’s electric field, not only in the first, but also in the second hydration shell. The hydration enthalpy was determined from the numerical simulation, taking into account corrections to the interfacial potential and to the spurious effects due to the periodicity imposed by the Ewald sums; the resulting value lied within the range of the various different experimental data. An analysis of the interaction energies between the ion and the water molecules in the different hydration shells and the bulk showed the same partition of the hydration enthalpy as for K⁺. The reason for this similarity is that at distances longer than 3 Å, the ion–water interaction is dominated by the charge-(enhanced) dipole term. Thus, it was concluded that starting at K⁺, the hydration properties of the heavier alkali metal cations should be very similar.     

Monte Carlo simulation of electrolytes in the constant voltage ensemble

The authors studied the structural, electrostatic, and electromechanical properties of the terlamellar structure composed of the anode, the cathode, and the electrolyte layer separating them. They used the Monte Carlo simulation technique in the constant voltage ensemble, where the electrical potential difference between the anode and the cathode is introduced as an external field. For ions, they used the primitive models of different sizes and valences in order to investigate how they affect the physical properties when an electrical field is applied between the electrodes. For electrodes, they used impermeable and permeable models, which mimic planar and porous electrodes, respectively. The asymmetry between the anions and the cations in size or valence was found to be responsible for the asymmetry in the concentration profile, the potential drop, and the stress distribution, in comparing the anode and the cathode sides. The charging/discharging process in the planar and porous electrodes is discussed at molecular level.     

Reactive Monte Carlo and grand-canonical Monte Carlo simulations of the propene metathesis reaction system

The influence of silicalite-1 pores on the reaction equilibria and the selectivity of the propene metathesis reaction system in the temperature range between 300 and 600 K and the pressure range from 0.5 to 7 bars has been investigated with molecular simulations. The reactive Monte Carlo (RxMC) technique was applied for bulk-phase simulations in the isobaric-isothermal ensemble and for two phase systems in the Gibbs ensemble. Additionally, Monte Carlo simulations in the grand-canonical ensemble (GCMC) have been carried out with and without using the RxMC technique. The various simulation procedures were combined with the configurational-bias Monte Carlo approach. It was found that the GCMC simulations are superior to the Gibbs ensemble simulations for reactions where the bulk-phase equilibrium can be calculated in advance and does not have to be simulated simultaneously with the molecules inside the pore. The confined environment can increase the conversion significantly. A large change in selectivity between the bulk phase and the pore phase is observed. Pressure and temperature have strong influences on both conversion and selectivity. At low pressure and temperature both conversion and selectivity have the highest values. The effect of confinement decreases as the temperature increases.