Although aqueous electrolytes are among the most important solutions, the molecular simulation of their intertwined properties of chemical potentials, solubility and activity coefficients has remained a challenging problem, and has attracted considerable recent interest. In this perspectives review, we focus on the simplest case of aqueous sodium chloride at ambient conditions and discuss the two main factors that have impeded progress. The first is lack of consensus with respect to the appropriate methodology for force field (FF) development. We examine how most commonly used FFs have been developed, and emphasise the importance of distinguishing between ‘Training Set Properties’ used to fit the FF parameters, and ‘Test Set Properties’, which are pure predictions of additional properties. The second is disagreement among solubility results obtained, even using identical FFs and thermodynamic conditions. Solubility calculations have been approached using both thermodynamic-based methods and direct molecular dynamics-based methods implementing coexisting solution and solid phases. Although convergence has been very recently achieved among results based on the former approach, there is as yet no general agreement with simulation results based on the latter methodology. We also propose a new method to directly calculate the electrolyte standard chemical potential in the Henry law ideality model. We conclude by making recommendations for calculating solubility, chemical potentials and activity coefficients, and outline a potential path for future progress.

We present an overview of novel numerical methods for chemical equilibrium and kinetic calculations for complex non-ideal multiphase systems. The methods we present for equilibrium calculations are based either on Gibbs energy minimization (GEM) calculations or on solving the system of extended law of mass-action (xLMA) equations. In both methods, no a posteriori phase stability tests, and thus no tentative addition or removal of phases during or at the end of the calculations, are necessary. All potentially stable phases are considered from the beginning of the calculation, and stability indices are immediately available at the end of the computation to determine which phases are actually stable at equilibrium. Both GEM and xLMA equilibrium methods are tailored for computationally demanding applications that require many rapid local equilibrium calculations, such as reactive transport modeling. The numerical method for chemical kinetic calculations we present supports both closed and open systems, and it considers a partial equilibrium simplification for fast reactions. The method employs an implicit integration scheme that improves stability and speed when solving the often stiff differential equations in kinetic calculations. As such, it requires compositional derivatives of the reaction rates to assemble the Jacobian matrix of the resultant implicit algebraic equations that are solved at every time step. We present a detailed procedure to calculate these derivatives, and we show how the partial equilibrium assumption affects their computation. These numerical methods have been implemented in Reaktoro (reaktoro.org), an open-source software for modeling chemically reactive systems. We finish with a discussion on the comparison of these methods with others in the literature.

This work reexamines seeded simulation results for NaCl nucleation from a supersaturated aqueous solution at 298.15 K and 1 bar pressure. We present a linear regression approach for analyzing seeded simulation data that provides both nucleation rates and uncertainty estimates. Our results show that rates obtained from seeded simulations rely critically on a precise driving force for the model system. The driving force vs. solute concentration curve need not exactly reproduce that of the real system, but it should accurately describe the thermodynamic properties of the model system. We also show that rate estimates depend strongly on the nucleus size metric. We show that the rate estimates systematically increase as more stringent local order parameters are used to count members of a cluster and provide tentative suggestions for appropriate clustering criteria.

The structures adopted within adsorbed monolayers of 4-n-octyl-4′-cyanobiphenyl (8-CB) molecules have been investigated using energy minimizations and molecular dynamics simulations of periodic systems. Using a smooth substrate potential, the most favourable energy of adsorption is found for a system with an eight-molecule unit-cell structure. This result is entirely consistent with scanning tunnelling microscopy studies of such systems, and differs from previous results using simulations of short strips which suggested a four-molecule unit cell. Molecular dynamics simulations of this 8-CB monolayer show that while the system exhibits smectic ordering at 150 K, the detailed eight-molecule unit-cell structure is lost. Simulations performed on a bilayer system indicate that the presence of a second molecular layer stabilizes the unit cell structure, except in the regions where there is partial penetration by the second layer molecules into the first layer. A third set of molecular dynamics simulations where the monolayer is confined between the substrate and a planar probe, shows that the eight-molecule unit cell is stable when out-of-plane motion is restricted by the probe. The effect of the molecular chain length on the intramolecular structure is also investigated: energy minimizations performed using the longer molecule 10-CB indicate that the eight-molecule unit cell is not the most stable configuration for this molecule. In this system, six- and ten-molecule unit cells both give lower energy arrangements than the eight-molecule cell adopted by 8-CB.

We present a molecular-level simulation study of CaCl₂ in water and crystalline hydrates formed by CaCl₂ 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. Moučka 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, CaCl₂·6H₂O) and at high-temperature conditions (sinjarite, CaCl₂·2H₂O). We implement our methodologies using for both phases the CaCl₂ 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⁻¹H₂O 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.