Novel barostat implementation for molecular dynamics.

We propose a novel implementation of the extended-dynamics equations for isothermal-isobaric ensemble in molecular dynamics, as the Martyna-Tobias-Klein thermostat and barostat. This method is suitable for systems with constraints and the Verlet-family integrators. Instead of iterations or the Trotter-expansion-based methods, both velocities and box sizes (scaling of bond lengths) are predicted. The algorithm begins with force calculation, requiring neither quarter nor half time steps, and necessitating iterations only inside SHAKE. Several tests demonstrate that the quality is comparable to other implementations. It is found that the formula relating the extended barostat mass to the characteristic time of volume fluctuations is inaccurate for condensed systems, which has consequences for the parameter setup. Emphasis is also put on the verification of the precise isothermal-isobaric ensemble and finite-size effects.

Properties of water and argon clusters developed in supersonic expansions.

Using adiabatic molecular dynamics coupled with the fluid dynamics equations, we model nucleation in an expanding beam of water vapor and argon on a microsecond scale. The size distribution of clusters, their temperature, and pickup cross sections in dependence on velocity are investigated and compared to the geometric cross sections and the experiment. The clusters are warmer than the expanding gas because of the time scale of relaxation processes. We also suggest that their translational and rotational kinetic energies are modified due to evaporative cooling. The pickup cross sections determined for the final clusters using molecules of the same kind increase with decreasing velocity, still obeying the (a+bN1/3)2 law.

Electron-induced fragmentation of water droplets: Simulation study.

The transport of free electrons in a water environment is still poorly understood. We show that additional insight can be brought about by investigating fragmentation patterns of finite-size particles upon electron impact ionization. We have developed a composite protocol aiming to simulate fragmentation of water clusters by electrons with kinetic energies in the range of up to 100 eV. The ionization events for atomistically described molecular clusters are identified by a kinetic Monte Carlo procedure. We subsequently model the fragmentation with classical molecular dynamics simulations, calibrated by non-adiabatic quantum mechanics/molecular mechanics simulations of the ionization process. We consider one-electron ionizations, energy transfer via electronic excitation events, elastic scattering, and also the autoionization events through intermolecular Coulombic decay. The simulations reveal that larger water clusters are often ionized repeatedly, which is the cause of substantial fragmentation. After losing most of its energy, low-energy electrons further contribute to fragmentation by electronic excitations. The simultaneous measurement of cluster size distribution before and after the ionization represents a sensitive measure of the energy transferred into the system by an incident electron.

Chemical potentials of alkaline earth metal halide aqueous electrolytes and solubility of their hydrates by molecular simulation: Application to CaCl2, antarcticite, and sinjarite.

We present a molecular-level simulation study of CaCl2 in water and crystalline hydrates formed by CaCl2 at ambient (298.15 K, 1 bar) conditions and at a high-temperature high-pressure state (365 K, 275 bars) typical of hydraulic fracturing conditions in natural-gas extraction, at which experimental properties are poorly characterized. We focus on simulations of chemical potentials in both solution and crystalline phases and on the salt solubility, the first time to our knowledge that such properties have been investigated by molecular simulation for divalent aqueous electrolytes. We first extend our osmotic ensemble Monte Carlo simulation technique [F. Moucka et al., J. Phys. Chem. B 115, 7849-7861 (2011)] to such solutions. We then describe and apply new methodology for the simulation of the chemical potentials of the experimentally observed crystalline hydrates at ambient conditions (antarcticite, CaCl2·6H2O) and at high-temperature conditions (sinjarite, CaCl2·2H2O). We implement our methodologies using for both phases the CaCl2 transferable force field (FF) based on simple point charge-extended water developed by Mamatkulov et al. [J. Chem. Phys. 138, 024505 (2013)], based on training sets involving single-ion and ion-pair low-concentration solution properties at near-ambient conditions. We find that simulations of the solution chemical potentials at high concentrations are somewhat problematic, exhibiting densities diverging from experimental values and accompanied by dramatically decreasing particle mobility. For the solid phases, the sinjarite crystalline lattice differs from experiment only slightly, whereas the simulations of antarcticite completely fail, due to instability of the crystalline lattice. The FF thus only successfully yields the sinjarite solubility, but its value m = 8.0(7) mol kg-1H2O lies well below the experimentally observed solubility range at 1 bar pressure of (12m, 15m) in the temperature interval (320 K, 400 K). We conclude that the used FF does not provide a good description of the experimental properties considered and suggest that improvement must take into account the crystalline properties.

Residual entropy of ice III from Monte Carlo simulation.

We calculated the residual entropy of ice III as a function of the occupation probabilities of hydrogen positions α and ß assuming equal energies of all configurations. To do this, a discrete ice model with Bjerrum defect energy penalty and harmonic terms to constrain the occupation probabilities was simulated by the Metropolis Monte Carlo method for a range of temperatures and sizes followed by thermodynamic integration and extrapolation to N = ∞. Similarly as for other ices, the residual entropies are slightly higher than the mean-field (no-loop) approximation. However, the corrections caused by fluctuation of energies of ice samples calculated using molecular models of water are too large for accurate determination of the chemical potential and phase equilibria.

Solubility of NaCl in water and its melting point by molecular dynamics in the slab geometry and a new BK3-compatible force field.

Saturated concentration of rock salt in water is determined by a simulation of brine in contact with a crystal in the slab geometry. The NaCl crystals are rotated to expose facets with higher Miller indices than [001] to brine. The rock salt melting point is obtained by both the standard and adiabatic simulations in the slab geometry with attention paid to finite size effects as well as to a possible influence of facets with higher Miller indices and applied stress. Two force fields are used, the Lennard-Jones-based model by Young and Cheatham with SPC/E water and the Kiss and Baranyai polarizable model with BK3 water. The latter model is refitted to thermomechanical properties of crystal NaCl leading to better values of solubility and the melting point.

Residual entropy of ices and clathrates from Monte Carlo simulation.

We calculated the residual entropy of ices (Ih, Ic, III, V, VI) and clathrates (I, II, H), assuming the same energy of all configurations satisfying the Bernal-Fowler ice rules. The Metropolis Monte Carlo simulations in the range of temperatures from infinity to a size-dependent threshold were followed by the thermodynamic integration. Convergence of the simulation and the finite-size effects were analyzed using the quasichemical approximation and the Debye-Hückel theory applied to the Bjerrum defects. The leading finite-size error terms, ln N/N, 1/N, and for the two-dimensional square ice model also 1/N(3/2), were used for an extrapolation to the thermodynamic limit. Finally, we discuss the influence of unequal energies of proton configurations.

Does the Sign of Charge Affect the Surface Affinity of Simple Ions?

The role the charge sign of simple ions plays in determining their surface affinity in aqueous solutions is investigated by computer simulation methods. For this purpose, the free surface of aqueous solutions of fictitious salts is simulated at finite concentration both with nonpolarizable point-charge and polarizable Gaussian-charge potential models. The salts consist of monovalent cations and anions that are, apart from the sign of their charge, identical to each other. In particular, we consider the small Na+ and the large I- ions together with their charge-inverted counterparts. In an attempt to avoid the interference even between the behavior of cations and anions, we also simulate systems containing only one of the above ions, and determine the free energy profile of these ions across the liquid-vapor interface of water at infinite dilution by potential of mean force (PMF) calculations. The obtained results reveal that, in the case of small ions, the anion is hydrated considerably stronger than the cation due to the close approach of water H atoms, bearing a positive fractional charge. As a consequence, the surface affinity of a small anion is even smaller than that of its cationic counterpart. However, considering that small ions are effectively repelled from the water surface, the importance of this difference is negligible. Further, a change in the hydration energy trends of the two oppositely charged ions is observed with their increasing size. This change is largely attributed to the fact that, with increasing ion size, the factor of 2 in the magnitude of the fractional charge of the closely approaching water atoms (i.e., O around cations and H around anions) outweighs the closer approach of the H than the O atom in the hydration energy. Thus, for large ions, being already surface active themselves, the surface affinity of the anion is larger than that of its positively charged counterpart. Further, such a difference is seen even in the case when the sign of the surface potential favors the adsorption of cations.

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.

Single-Particle Dynamics at the Intrinsic Surface of Aqueous Alkali Halide Solutions.

The distribution of ions in the proximity of the liquid-vapor interface of their aqueous solution has been the subject of an intense debate during the last decade. The effects of ionic polarizability have been one of its salient aspects. Much less has been said about the corresponding dynamical properties, which are substantially unexplored. Here, we investigate the single-particle dynamics at the liquid-vapor interface of several alkali halide solutions, using molecular dynamics simulations with polarizable and nonpolarizable force fields and intrinsic surface analysis. We analyze the diffusion coefficient, residence time, and velocity autocorrelation function of water and ions and investigate how these properties depend on the molecular layer where they reside. While anions are found in the first molecular layer for relatively long times, cations are only making quick excursions into it, thanks to thermal fluctuations. The in-layer residence time of ions and their molar fraction in the layer turned out to be linearly dependent on each other. We interpret this unexpected result using a simple two-state model. In addition, we found that, unlike water and other neat molecular liquids that show a different diffusion mechanism at the surface than in the bulk of their liquid phase, ions do not enjoy enhanced mobility in the surface layer of their aqueous solution. This result indicates that ions in the surface layer are shielded by their nearest water neighbors from being exposed to the vapor phase as much as possible. Such positions are available for the ions at the negatively curved troughs of the molecularly rugged liquid surface.

Molecular Dynamics of Heterogeneous Systems on GPUs and Their Application to Nucleation in Gas Expanding to a Vacuum.

Expansion of water vapor through a small orifice to a vacuum produces liquid or frozen clusters which in the experiment serve as model particles for atmospheric aerosols. Yet, there are controversies about the shape of these clusters, suggesting that the nucleation process is not fully understood. Such questions can be answered by molecular dynamics simulations; however, they require microsecond-scale runs with thousands of molecules and accurate energy conservation. The available highly parallel codes typically utilize domain decomposition and are inefficient for heterogeneous systems as clusters in a dilute gas. In this work, we present an implementation of molecular dynamics on graphics processing units based on the Verlet list and apply it to several systems for which experimental data are available. We reproduce sufficiently sized clusters but not the experimentally observed clusters of irregular shape.

Virial coefficients and equation of state of the penetrable sphere model.

We study the penetrable sphere (alias square mound) model in the fluid phase by means of the virial expansion, molecular dynamics simulations, and Ornstein-Zernike integral equation. The virial coefficients up to B(8) are expressed as polynomials in the Boltzmann factor with the coefficients calculated by a Monte Carlo integration. New data for pressure and internal energy are obtained by molecular dynamics simulations with attention paid to finite-size errors and properties of the Andersen thermostat. The data and virial coefficients are correlated by a formula for the Helmholtz free energy. We also propose a new closure for the Ornstein-Zernike equation and test several other closures.

Pressure in Molecular Simulations with Scaled Charges. 1. Ionic Systems.

Charge scaling, rationalized as MDEC (molecular dynamics in electronic continuum) or ECC (electronic continuum correction), has become a widely used simple approach to how to avoid self-consistent induced dipoles yet approximately take into account the effects of electronic polarizability. It has been assumed that the continuum permittivity does not depend on density; in turn, pressure is calculated by standard formulas. In this work, we elaborate a complementary approximation of density-independent molecular polarizability and derive formulas for pressure corrections within the MDEC framework; real behavior lies between these two extremes. The pressure corrections for test ionic systems are huge and negative, leading to sizable densities in constant-pressure MDEC simulations. A comparison of MDEC results with equivalent polarizable systems gives a good pressure match for a crystal but very low MDEC pressures for ionic liquids. These results witness about the importance of a correct density dependence not only of continuum permittivity in MDEC simulations but also of polarizability in polarizable simulations.

Surface Affinity of Alkali and Halide Ions in Their Aqueous Solution: Insight from Intrinsic Density Analysis.

The surface tension of all aqueous alkali halide solutions is higher than that of pure water. According to the Gibbs adsorption equation, this indicates a net depletion of these ions in the interfacial region. However, simulations and experiments show that large, soft ions, such as I-, can accumulate at the liquid/vapor interface. The presence of a loose hydration shell is usually considered to be the reason for this behavior. In this work, we perform computer simulations to characterize the liquid-vapor interface of aqueous alkali chloride and sodium halide solutions systematically, considering all ions from Li+ to Cs+ and from F- to I-. Using computational methods for the removal of surface fluctuations, we analyze the structure of the interface at a dramatically enhanced resolution, showing that the positive excess originates in the very first molecular layer and that the next 3-4 layers account for the net negative excess. With the help of a fictitious system with charge-inverted ion pairs, we also show that it is not possible to rationalize the surface affinity of ions in solutions in terms of the properties of anions and cations separately. Moreover, the surface excess is generally dominated by the smaller of the two ions.

Free Energy of Classical Molecular Crystals by Thermodynamic Integration from a Harmonic Reference.

We develop an algorithm for calculating the normal modes of vibration of mechanical systems with constraints, particularly of molecules with rigid bonds and models of rigid molecules, and use it to obtain the harmonic free energy of a crystal. The anharmonic correction is then calculated by the conventional thermodynamic integration over temperature in the NVT ensemble. Attention is paid to finite-size errors, tail corrections, thermostat choice, ergodicity, and other sources of inaccuracies. The calculated free energy of ice XIV modeled by the TIP4P/2005 potential agrees with the previously reported value and is by one order more accurate.

Direct Molecular Dynamics Simulation of Nucleation during Supersonic Expansion of Gas to a Vacuum.

We develop a methodology for direct molecular-level simulation of adiabatic expansion of gas through a small orifice to a vacuum. The gas attains supersonic speeds, cools, and nucleates. The proposed approach combines equations of frictionless fluid dynamics with molecular dynamics simulation in an expanding periodic box. There are two key components of the proposed algorithm: (i) a time-reversible integrator tailored to an expanding system, and (ii) an iterative procedure employed to satisfy the condition of steady flow. For a conical nozzle (opening angle of 60°), the simulations with argon and water vapor predict cluster sizes in agreement with the experiment. Clusters of irregular shapes observed in the experiment [J. Lengyel et al. Phys. Rev. Lett. 2014, 112, 113401] are not reproduced. The role of friction, turbulence, and sonic boom originating at the sharp nozzle edge is discussed.

Accurate Binding of Sodium and Calcium to a POPC Bilayer by Effective Inclusion of Electronic Polarization.

Binding affinities and stoichiometries of Na+ and Ca2+ ions to phospholipid bilayers are of paramount significance in the properties and functionality of cellular membranes. Current estimates of binding affinities and stoichiometries of cations are, however, inconsistent due to limitations in the available experimental and computational methods. In this work, we improve the description of the binding details of Na+ and Ca2+ ions to a 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) bilayer by implicitly including electronic polarization as a mean field correction, known as the electronic continuum correction (ECC). This is applied by scaling the partial charges of a selected state-of-the-art POPC lipid model for molecular dynamics simulations. Our improved ECC-POPC model reproduces not only the experimentally measured structural parameters for the ion-free membrane, but also the response of lipid headgroup to a strongly bound cationic amphiphile, as well as the binding affinities of Na+ and Ca2+ ions. With our new model, we observe on the one side negligible binding of Na+ ions to POPC bilayer, while on the other side stronger interactions of Ca2+ primarily with phosphate oxygens, which is in agreement with the previous interpretations of the experimental spectroscopic data. The present model results in Ca2+ ions forming complexes with one to three POPC molecules with almost equal probabilities, suggesting more complex binding stoichiometries than those from simple models used to interpret the NMR data previously. The results of this work pave the way to quantitative molecular simulations with realistic electrostatic interactions of complex biochemical systems at cellular membranes.

Molecular simulations of aqueous electrolyte solubility: 1. The expanded-ensemble osmotic molecular dynamics method for the solution phase.

We have developed a molecular-level simulation technique called the expanded-ensemble osmotic molecular dynamics (EEOMD) method, for studying electrolyte solution systems. The EEOMD method performs simulations at a fixed number of solvent molecules, pressure, temperature, and overall electrolyte chemical potential. The method combines elements of constant pressure-constant temperature molecular dynamics and expanded-ensemble grand canonical Monte Carlo. The simulated electrolyte solution systems contain, in addition to solvent molecules, full and fractional ions and undissociated electrolyte molecular units. The fractional particles are coupled to the system via a coupling parameter that varies between 0 (no interaction between the fractional particle and the other particles in the system) and 1 (full interaction between the fractional particle and the other particles in the system). The time evolution of the system is governed by the constant pressure-constant temperature equations of motion and accompanied by random changes in the coupling parameter. The coupling-parameter changes are accepted with a probability derived from the expanded-ensemble osmotic partition function corresponding to the prescribed electrolyte chemical potential. The coupling-parameter changes mimic insertion/deletion of particles as in a crude grand canonical Monte Carlo simulation; if the coupling parameter becomes 0, the fractional particles disappear from the system, and as the coupling parameter reaches unity, the fractional particles become full particles. The method is demonstrated for a model of NaCl in water at ambient conditions. To test our approach, we first determine the chemical potential of NaCl in water by the thermodynamic integration technique and by the expanded-ensemble method. Then, we carry out EEOMD simulations for different specified values of the overall NaCl chemical potential and measure the concentration of ions resulting from the simulations. Both computations give consistent results, validating the EEOMD methodology.

Virial coefficients of hard spheres and hard disks up to the ninth.

A technique for topological analysis of the Ree-Hoover diagrams is developed with the aim to calculate the Ree-Hoover weights up to the ninth order with moderate demands on computer storage and CPU time. The ninth virial coefficients of hard spheres and disks are calculated, and the lower virial coefficients are accurately recalculated. The calculations require several spanning diagrams; the most important spanning chains are generated by reptation, other spanning diagrams by the standard Metropolis Monte Carlo algorithm. The tenth and eleventh virial coefficients for hard spheres are estimated.

Static Dielectric Constant from Simulations Revisited: Fluctuations or External Field?

We compare two methods-the fluctuation formula and the application of an external electric field-for determining the static relative permittivity (static dielectric constant) from molecular simulations in the Ewald dielectric boundary conditions. We express both the systematic and statistical errors in terms of a dimensionless saturation (dielectric polarization reduced by its maximum) and show that both methods possess the same efficiency. Saturation in the fluctuation route depends on the number of particles and the permittivity of the surrounding medium, where a value of infinity (tinfoil) for the latter is usually reasonable but not necessarily optimum. Saturation in the external field route is easily controlled by the field intensity; hence, we develop an extrapolation scheme and apply it to several water models, specifically, nonpolarizable (SPC/E, TIP4P/2005, NE6) and polarizable (POL3, DC97, BK3). In addition, we discuss the Clausius-Mossotti approximation for the optical permittivity; this quantity is needed in the fluctuation formula for polarizable models.