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.     

A classical polarizable model for simulations of water and ice

We develop a classical rigid polarizable model of water for molecular simulations of water and ice. The model uses the Rowlinson five-site geometry: oxygen bearing the Lennard-Jones interaction and linearly polarizable point dipole, two positively charged hydrogens, and two massless negative charges placed symmetrically off oxygen so that the experimental dipole moment is reproduced. The target properties are the densities of water and ice Ih, diffusivity, enthalpies of fusion and vaporization, and the ice Ih melting point. The surface tension at lower temperatures is by 7% underestimated whereas the dielectric constant by 6% overestimated. Diffusivity and viscosity worsen at higher temperatures, although the Stokes radius is overestimated only by 2-7%. The ice Ih melting temperature is 260 K and the temperature of maximum density is 269 K. Rescaling the charges by a factor of 1.01 and Lennard-Jones energy by 1.0201 improves the melting point and energy-related quantities but shifts the agreement of kinetic properties to higher temperatures. For the model we propose abbreviation POL4D.     

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.     

An efficient parallelization scheme for molecular dynamics simulations with many-body, flexible, polarizable empirical potentials: application to water

An efficient parallelization scheme for classical molecular dynamics simulations with flexible, polarizable empirical potentials is presented. It is based on the standard Ewald summation technique to handle the long-range electrostatic and induction interactions. The algorithm for this parallelization scheme is designed for systems containing several thousands of polarizable sites in the simulation box. Its performance is evaluated during molecular dynamics simulations under periodic boundary conditions with unit cell sizes ranging from 128 to 512 molecules employing two flexible polarizable water models [DC(F) and TTM2.1-F] containing 1 and 3 polarizable sites, respectively. The time-to-solution for these two polarizable models is compared with the one for a flexible, pairwise-additive water model (TIP4F). The benchmarks were performed on both shared and distributed memory platforms. As a result of the efficient calculation of the induced dipole moments, a superlinear scaling as a function of the number of the processors is observed. To the best of our knowledge, this is the first reported results of parallel scaling and performance for simulations of liquid water with a polarizable potential under periodic boundary conditions.     

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.     

Monte Carlo simulation of polarizable systems: Early rejection scheme for improving the performance of adiabatic nuclear and electronic sampling Monte Carlo simulations

An early rejection scheme for trial moves in adiabatic nuclear and electronic sampling Monte Carlo simulation (ANES-MC) of polarizable intermolecular potential models is presented. The proposed algorithm is based on Swendsen–Wang filter functions for prediction of success or failure of trial moves in Monte Carlo simulations. The goal was to reduce the amount of calculations involved in ANES-MC electronic moves, by foreseeing the success of an attempt before making those moves. The new method was employed in Gibbs ensemble Monte Carlo (GEMC) simulations of the polarizable simple point charge-fluctuating charge (SPC-FQ) model of water. The overall improvement in GEMC depends on the number of swap attempts (transfer molecules between phases) in one Monte Carlo cycle. The proposed method allows this number to increase, enhancing the chemical potential equalization. For a system with 300 SPC-FQ water molecules, for example, the fractions of early rejected transfers were about 0.9998 and 0.9994 at 373 and 423 K, respectively. This means that the transfer moves consume only a very small part of the overall computing effort, making GEMC almost equivalent to a simulation in the canonical ensemble.     

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     

Chain conformation and solvent partitioning in reversed-phase liquid chromatography: Monte Carlo simulations for various water/methanol concentrations

Many structural models for the stationary phase in reversed-phase liquid chromatography (RPLC) systems have been suggested from thermodynamic and spectroscopic measurements and theoretical considerations. To provide a molecular picture of chain conformation and solvent partitioning in a typical RPLC system, a particle-based Monte Carlo simulation study is undertaken for a dimethyl octadecyl (C(18)) bonded stationary phase on a model siliceous substrate in contact with mobile phases having different methanol/water concentrations. Following upon previous simulations for gas-liquid chromatography and liquid-liquid phase equilibria, the simulations are conducted using the configurational-bias Monte Carlo method in the Gibbs ensemble and the transferable potentials for phase equilibria force field. The simulations are performed for a chain surface density of 2.9 micromol/m(2), which is a typical bonded-phase coverage for mono-functional alkyl silanes. The solvent concentrations used here are pure water, approximately 33 and 67% mole fraction of methanol and pure methanol. The simulations show that the chain conformation depends only weakly on the solvent composition. Most chains are conformationally disordered and tilt away from the substrate normal. The interfacial width increases with increasing methanol content and, for mixtures, the solvent shows an enhancement of the methanol concentration in a 10 Angstrom region outside the Gibbs dividing surface. Residual surface silanol groups are found to provide hydrogen bonding sites that lead to the formation of substrate bound water and methanol clusters, including bridging clusters that penetrate from the solvent/chain interfacial region all the way to the silica surface.     

Monte Carlo simulations of equilibrium solubilities and structure of water in n-alkanes and polyethylene

Gibbs ensemble Monte Carlo methods based on a force field that combines the simple point charge [Berendsen et al., in Intermolecular Forces, edited by Pullman (Reidel, Dordrecht, 1981), p. 331] and transferable potentials for phase equilibria [Martin and Siepmann, J. Phys. Chem. B 102, 2659 (1998)] models were used to study the equilibrium properties of binary systems consisting of water and n-alkanes with chain lengths from hexane to hexadecane. In addition, systems where extended linear alkane chains (up to 300 carbon units long) were used to represent amorphous polyethylene were simulated in the presence of water using a connectivity altering osmotic Gibbs ensemble. In these simulations the equilibrium between a liquid water phase and a polymer phase into which water was inserted was studied. The predicted solubilities, which were determined between 350 and 550 K, are in good agreement with experiment, where experimental results are available, and the density of water molecules in the hydrocarbons is approximately 63% as high as in saturated water vapor under the same conditions. At the lower temperatures most of the water exists as monomers; increasing the temperature leads to an increase in the density of water in the alkane phase and hence in the fraction of molecules that participate in clusters. Dimers are the most prevalent clusters in all hydrocarbons and at all temperatures studied, and the fraction of clusters of given size decrease with increasing cluster size. A large fraction of trimers, tetramers, and pentamers, which are the cluster sizes for which topologies have been studied, are cyclic at low temperatures, but at higher temperatures linear structures predominate. The same properties are observed for pure water vapor clusters in equilibrium with the liquid phase, showing that the cluster topologies are not significantly affected by the surrounding hydrocarbon.     

Macroion solutions in the cell model studied by field theory and Monte Carlo simulations

Aqueous solutions of charged spherical macroions with variable dielectric permittivity and their associated counterions are examined within the cell model using a field theory and Monte Carlo simulations. The field theory is based on separation of fields into short- and long-wavelength terms, which are subjected to different statistical-mechanical treatments. The simulations were performed by using a new, accurate, and fast algorithm for numerical evaluation of the electrostatic polarization interaction. The field theory provides counterion distributions outside a macroion in good agreement with the simulation results over the full range from weak to strong electrostatic coupling. A low-dielectric macroion leads to a displacement of the counterions away from the macroion.     

Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice I(h)

We present a molecular dynamics simulation study in which we determined the melting point of ice I(h) for the polarizable SWM4-NDP water model (Lamoureux et al., Chem. Phys. Lett., 2006, 418, 245-249) and compared the performance of several popular water force fields, both polarizable and nonpolarizable, in terms of melting temperature, stability and orientational structuring of ice. The simulations yield the melting temperature of SWM4-NDP ice as low as T(m) = 185 ± 10 K, despite the quadrupole moment of a SWM4-NDP water molecule being close to the experimental gas phase value. The results thus show that the dependence of T(m) on the molecular quadrupole, observed for the three- and four-site water models, is generally lost if polarization is explicitly included. The study also shows that adding polarizability to a planar three-charge water model increases orientational disorder in hexagonal ice. In addition, analysis of the tetrahedral order in bulk ice reveals a correlation between the pre-existing degree of orientational disorder in ice simulated using different polarizable and nonpolarizable models and the melting temperature of the models. Our findings thus suggest some new considerations regarding the role of polarization forces in a crystalline solid that may guide future development of reliable polarizable water models for ice.     

Monte Carlo simulations of gas solubility in the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate

The Henry's constants of water, carbon dioxide, ethane, ethene, methane, oxygen, and nitrogen are computed in the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate ([bmim][PF(6)]) using test particle insertion and expanded ensemble Monte Carlo methods. The partial molar enthalpy and partial molar entropy of solvation are also computed for water, carbon dioxide, and oxygen. The results from the simulations are compared against experimental data from the literature. In addition, the accuracy and precision of the two methods in determining the Henry's constant are examined. Local organization of the ionic liquid around a solute molecule is analyzed, and the interactions responsible for the experimentally observed solubility trends are identified.     

Monte Carlo simulations of solid nitrogen in the isothermal–isobaric ensemble with an ab initioSCF–CI potential

Monte Carlo simulations in the isothermal–isobaric ensemble for the α phase of solid N₂ have been carried out with two different pair potentials obtained from ab initio quantum chemical calculations. Comparison is made with data obtained from empirical potentials as well as with experimental data, and in general, the agreement with experiment is found to be good. It is also found that the differences between experiment and theory can largely be attributed to correlation effects, not considered in the quantum mechanical calculations of the pair potential.     

Thermodynamic properties of model solids with short-ranged potentials from Monte Carlo simulations and perturbation theory

Monte Carlo simulations have been performed to determine the excess energy and the equation of state of fcc solids with Sutherland potentials for wide ranges of temperatures, densities, and effective potential ranges. The same quantities have been determined within a perturbative scheme by means of two procedures: (i) Monte Carlo simulations performed on the reference hard-sphere system and (ii) second-order Barker-Henderson perturbation theory. The aim was twofold: on the one hand, to test the capability of the "exact" MC-perturbation theory of reproducing the direct MC simulations and, on the other hand, the reliability of the Barker-Henderson perturbation theory, as compared with direct MC simulations and MC-perturbation theory, to determine the thermodynamic properties of these solids depending on temperature, density, and potential range. We have found that the simulation data for the excess energy obtained from the two procedures are in close agreement with each other. For the equation of state, the results from the MC-perturbation procedure also agree well with the direct MC simulations except for very low temperatures and extremely short-ranged potentials. Regarding the Barker-Henderson perturbation theory, we have found that in general the second-order approximation does not provide significant improvement over the first-order one.     

Attraction of iodide ions by the free water surface, revealed by simulations with a polarizable force field based on Drude oscillators

Recent theoretical and experimental studies have shown that polarizable anions, such as iodide and bromide, preferentially accumulate close to the surface of electrolyte solutions. This finding is in sharp contrast to the previously prevailing idea that salts are dielectrically excluded from the free water surface and opens up new avenues for research in specific salt effects. In this work, we have verified the ability of a recently introduced polarizable water model, SWM4-DP, to reproduce this behavior, by simulations of a NaI/water slab, corresponding to a 1.2 M solution. The water and ion polarizabilities are modeled by classical Drude oscillator particles. As revealed by the simulations, a double layer is formed close to the free water surface, with the iodide ions located closer to the interface and the sodium ions at a neighboring, interior layer. Near the surface, all solution species acquire an induced dipole moment, that is perpendicular to the surface and points toward the exterior. The double charge layer causes ordering of water at a subsurface region. Simulations with a simpler system of a single iodide ion in a water slab show that the surface position is stabilized by induced charge interactions; in contrast, the charge-dipole interactions between the iodide permanent charge and the water permanent dipole moment favor the bulk position. Thus, the polarizabilities of ion and water are essential for explaining the increased preference of iodide for the air-water interface, in accordance with other studies.     

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.     

Phase separation in mixtures of Yukawa and charged Yukawa particles from Gibbs ensemble Monte Carlo simulations and the mean spherical approximation

The phase equilibrium of mixtures of Yukawa and charged Yukawa particles is studied by means of Gibbs ensemble Monte Carlo (GEMC) simulation method and the mean spherical approximation (MSA). The strength of the Coulomb energy compared to that of the Yukawa attraction is characterized by a coupling constant. For low coupling constants a classical vapor--liquid phase separation appears with a good agreement between GEMC and the MSA. For high coupling constant, a phase separation between a salt poor and a salt rich phase occurs that resembles the phase equilibrium behavior of the solvent primitive model.     

Can Monte Carlo detect the absence of ordering in a model liquid crystal?

We investigate the possibility of inferring the absence of an ordered phase using Monte Carlo simulations. The example we have chosen is that of a one dimensional Lebwohl-Lasher model, where an analytic solution is available. We argue that Monte Carlo can be of help even in this delicate sector notwithstanding the complications created by periodic boundary conditions.