This article reviews recent molecular simulation studies of “collective” properties of aqueous electrolyte solutions, specifically free energies and activity coefficients, solubilities, nucleation rates of crystals, and transport coefficients. These are important fundamental properties for biology and geoscience, but also relevant for many technological applications. Their determination from molecular-scale calculations requires large systems and long sampling times, as well as specialized sampling algorithms. As a result, such properties have not typically been taken into account during optimization of force field parameters; thus, they provide stringent tests for the transferability and range of applicability of proposed molecular models. There has been significant progress on simulation algorithms to enable the determination of these properties with good statistical uncertainties. Comparisons of simulation results to experimental data reveal deficiencies shared by many commonly used models. Moreover, there appear to exist specific tradeoffs within existing modeling frameworks so that good prediction of some properties is linked to poor prediction for specific other properties. For example, non-polarizable models that utilize full charges on the ions generally fail to predict accurately both activity coefficients and solubilities; the concentration dependence of viscosity and diffusivity for these models is also incorrect. Scaled-charge models improve the dynamic properties and could also perform well for solubilities but fail in the prediction of nucleation rates. Even models that do well at room temperature for some properties generally fail to capture their experimentally observed temperature dependence. The main conclusion from the present review is that qualitatively new physics will need to be incorporated in future models of electrolyte solutions to allow the description of collective properties for broad ranges of concentrations, temperatures, and solvent conditions.

Electrolyte solutions play important roles in many technological applications, e.g., separation and purification operations and energy production. They provide the medium in which biological functions take place and constitute a major component of the geological environment. The thermodynamic and transport properties of electrolyte solutions are strongly influenced by the long-range Coulombic interactions. Of particular importance for many applications are “collective” properties, such as activity coefficients, solubilities, nucleation rates of crystals from solutions, and transport coefficients, e.g., ionic conductivities. From a microscopic point of view, these are properties that require sampling of relatively large systems over time scales covering the relevant relaxation times, which typically are much longer than molecular vibration time scales. Sophisticated sampling algorithms are needed to obtain statistical quantities (such as the free energies) that underlie many of the collective properties. These properties generally depend on temperature, in many cases strongly so. In these respects, collective properties stand in contrast to simpler structural or thermodynamic properties, such as pair correlation functions or bulk densities, which have shorter relaxation times, a generally weaker dependence on temperature, and for which there are straightforward ways to obtain them directly from simulation trajectories that contain a limited number of independent molecular configurations. As a result of these difficulties, the collective properties of electrolyte solutions have not been extensively studied by simulations until relatively recently, and they are typically not included as target properties in the parameterization of water and electrolyte force fields.

A large number of careful experimental measurements of the properties of electrolyte solutions are available, primarily for aqueous solutions of single salts near room temperature. Thermodynamic non-idealities in electrolyte solutions can be described at low concentrations using the Debye–Hückel theory, but deviations become significant above 0.1M for 1:1 electrolytes, and at even lower concentrations for multivalent salts.1 Several phenomenological models have been proposed to describe activity coefficients in concentrated solutions of electrolytes,2 but these generally rely on a number of adjustable parameters fitted to experimental data. Thus, there is a real need for modeling methods that can be used to obtain collective properties over a broad range of compositions, temperatures, and pressures.

An emerging approach for describing the properties of electrolyte solutions involves the use of computational methods. Presently, quantum mechanical methods that take into account electronic and nuclear degrees of freedom are not able to cover the time and length scales needed for sampling collective properties. A more practical alternative is the use of molecular simulation methods utilizing classical molecular dynamics or Monte Carlo sampling. These methods require as input “force fields” consisting of functional forms (e.g., the Lennard-Jones 12-6 expression) and a specific choice of model parameters (e.g., well depth and size in the case of the Lennard-Jones expression) describing the interactions between molecules in a system.

The focus of the current article is methodological developments of recent years that enable the computation of collective properties with good statistical accuracy from classical molecular dynamics and Monte Carlo methods using explicit-solvent force fields. Implicit-solvent (or “primitive”) models that replace the solvent with a dielectric continuum were quite popular in early computer simulations of electrolyte solutions3 because of their simplicity and speed but cannot describe the hydration structure or changes in the properties of the solution at higher concentrations. We also do not consider here coarse-grained models4 often used for large-scale biomolecular simulations.

Because collective properties have not typically been used in developing force fields for electrolytes, they provide rather stringent tests for their transferability and range of applicability. Such considerations and tests are fortunately becoming more widely accepted.5 As was previously the case when methods became available to obtain collective properties for simple fluids and mixtures (e.g., vapor–liquid coexistence and critical parameters6,7), we find that calculations of collective properties made possible by algorithmic advances reveal previously unidentified weaknesses that need to be remedied in future generations of force fields. We restrict our attention here to fully dissociated simple ions, such as alkali and alkali earths or halides, in aqueous solution. In addition, in order to make comparisons across different properties and to keep the size of this article manageable, we focus primarily on systems for which two important collective properties, namely, the activity coefficients and solubilities, are available. There are many other systems or models for which only a limited subset of properties have been calculated, which will not be covered here. The importance of solution chemical potentials and solubilities in evaluating force fields for electrolytes has been clearly pointed out, e.g., by Aragones et al.8 and Moučka et al.9 

There have been several prior articles that overlap in scope with the present review. Joung and Cheatham tested several ion and water model combinations using activity coefficient and solubility calculations.10 Moučka et al. compared several common NaCl force fields with the SPC/E water model in terms of chemical potentials and solubilities.9 The thermodynamic and transport properties of aqueous NaCl solutions were discussed by Jiang et al.11 Two recent comprehensive reviews on explicit-solvent simulations of aqueous solutions for the calculation of chemical potentials, solubilities, and activity coefficients have been presented by Nezbeda and co-workers.12,13

The structure of this “perspective” article is as follows. The article starts with a brief review of molecular models for water and simple ions in solution relevant for the subsequent discussion of simulations of collective properties. Recent methodological advances in simulation methods are presented next, along with the results obtained using molecular-based, explicit-solvent models. The material is organized on the basis of different properties, starting with activity coefficients and solubilities, and followed by transport coefficients and nucleation rates. Comparisons to experimental data are made for all properties of interest. While many of the properties are described reasonably well, common patterns emerge on model deficiencies that need to be addressed by future force field development, incorporating additional features absent from the current generation of models.

A key consideration in explicit-solvent classical simulations of electrolyte solutions is the force field for water, which typically represents the majority of molecules in a simulated system. The current state-of-the-art water models contain many-body potential energy terms and can describe correctly the structure and dynamics of small water clusters, the liquid, and even ice phases.14 However, because of their computational cost, these sophisticated many-body models have not been used in combination with advanced sampling algorithms to determine collective properties such as the free energy or phase behavior. For example, for the MB-pol “first principles” water model of Babin and co-workers,15–17 even the relatively simple-to-compute vapor–liquid phase behavior of the pure fluid has not yet been determined.

The most widely used water force fields in simulations are fixed-point-charge models, in the same family as the early ST2 model of Stillinger and Rahman.18 The comprehensive review by Vega and Abascal19 examined several rigid, non-polarizable models in this family and found that the TIP4P/2005 model20 performed best on the basis of deviations from experimental data for a specific set of properties, primarily of the liquid. Interestingly, the TIP3P model21 popular for biomolecular simulations performed worse among the models studied. Shvab and Sadus22 examined a broader set of models and concluded that, near ambient conditions, ab initio inspired, many-parameter models such as iAMOEBA and MCYna performed the best for the properties examined, followed by polarizable models such as BK3.23 Unfortunately, none of the ab initio inspired models have been studied with respect to their ability to model the collective properties of electrolyte solutions to date.

The selected water models are summarized in Table I—the basis for the selection is the availability of both activity coefficients and solubilities for specific electrolyte models, as discussed in Secs. III–VI. This is, of course, a tiny fraction of the available models for water. This table shows the primary reference and year in which it appeared, the model name (or abbreviation used in the present article), types of non-Coulombic interactions used for the description of dispersion and repulsion interactions between the oxygen atoms on different molecules, the use of point charges or Gaussian charge distributions for the description of the Coulombic interactions, and treatment of polarizability, if any. All of the models listed in Table I have rigid geometries. The two options for the description of the non-Coulombic interactions are the Lennard-Jones (LJ) 12-6 and exponential-6 (Buckingham) functional forms. Coulombic interactions in the selected models are described either as fixed point charges at specified locations, or Gaussian charge distributions, as indicated in this table. Polarizability in the BK3 model is treated through the use of charge-on-spring Drude oscillators.30 Model parameters (such as potential well depths, size parameters, detailed geometries, and partial charges) can be obtained from the original references.

TABLE I.

Selected water models.

ReferencesYearNameNon-CoulombicCoulombicPolarizability
21  1983 TIP3P LJ Points  
24  1987 SPC/E LJ Points  
25  1998 Exp-6 Exponential-6 Points  
26  2003 SWM4 LJ Points Drude 
27  2004 TIP4Pew LJ Points  
20  2005 TIP4P/2005 LJ Points  
23  2013 BK3 Exponential-6 Gaussian Drude 
28  2014 TIP4P/ϵ LJ Points  
29  2015 SPC/ϵ LJ Points  
ReferencesYearNameNon-CoulombicCoulombicPolarizability
21  1983 TIP3P LJ Points  
24  1987 SPC/E LJ Points  
25  1998 Exp-6 Exponential-6 Points  
26  2003 SWM4 LJ Points Drude 
27  2004 TIP4Pew LJ Points  
20  2005 TIP4P/2005 LJ Points  
23  2013 BK3 Exponential-6 Gaussian Drude 
28  2014 TIP4P/ϵ LJ Points  
29  2015 SPC/ϵ LJ Points  

Most existing ion force fields were developed for use in combination with specific fixed-point-charge water models. There are relatively few ion force fields compatible with polarizable water models and almost none developed for use with ab initio inspired many-body potentials for water. Typical properties used in the development of ion models include hydration energies in solution and lattice energies and spacing of the ionic crystals. Some models were developed taking into account solution chemical potentials31 or the derivatives of the chemical potential with respect to concentration through the Kirkwood–Buff integral formalism.32 Monovalent alkali metal and halide ions, most prominently NaCl, have been the target of significant force field development and are the systems for which a broad set of results for the collective properties of interest in the present review are available. Multivalent electrolytes exhibit stronger deviations from ideality at lower concentrations. Because of the stronger interactions, they present additional sampling difficulties relative to the monovalent case. Their solid phases are frequently hydrated with varying ratios of water to ions stable at different temperatures.

The suggestion has been made33,34 that common ion models overestimate interactions in solution by considering ions as having full charges, as one would naturally expect for strong electrolytes that fully dissociate in solution. Using full charges effectively ignores electronic polarization effects of the solvent and is likely one of the reasons underlying the observation that practically all models for which salt solubilities in water have been determined underestimate the experimental values, as described later in this article. A proposed solution for this issue is to scale the charges by the inverse square root of the electronic part of the dielectric constant of water. Such charge scaling is routinely used in simulations of ionic liquids and in biophysical simulations. However, as will be discussed in Secs. III–VI of the present article, charge scaling entails compromises with respect to other collective properties, e.g., activity coefficients at low-to-intermediate concentrations and nucleation rates.

Table II summarizes the ion force fields for which data for the collective properties of interest in the context of the present article are directly available. Many, but not all, of the force fields listed in Table II are for NaCl, which has been the focus of the majority of studies of collective properties. This table lists the primary reference and year proposed, the model name, water models for which the ion force field has been optimized, and the types of interactions present in the model. Specifically, non-Coulombic interactions in the models of interest are described using the LJ, exponential-6 (Buckingham), or Born–Huggins–Mayer functional forms. The table indicates if models use point charges or Gaussian charge distributions for the Coulombic interactions, the inclusion of polarizability (if any), and the use of scaled or full charges.

TABLE II.

Selected ion models.

ReferencesYearNameH2O model(s)Non-CoulombicCoulombicPolarizabilityCharges
35  1964 TF  Born–Huggins–Mayer Points  Full 
36  1994 SD SPC/E LJ Points  Full 
37  2006 AH/SWM4 SWM4 LJ Points Drude Full 
38  2008 JC TIP3P, TIP4Pew, SPC/E LJ Points  Full 
32  2011 KBFF SPC/E LJ Points  Full 
39  2012 MP  LJ Points  Scaled 
40  2014 RDVH SPC/E LJ Points  Full 
41  2014 AH/BK3 BK3 Exponential-6 Gaussian Drude Full 
42  2016 NaCl/ϵ TIP4P/ϵ, SPC/ϵ LJ Points  Scaled 
43  2017 Madrid TIP4P/2005 LJ Points  Scaled 
ReferencesYearNameH2O model(s)Non-CoulombicCoulombicPolarizabilityCharges
35  1964 TF  Born–Huggins–Mayer Points  Full 
36  1994 SD SPC/E LJ Points  Full 
37  2006 AH/SWM4 SWM4 LJ Points Drude Full 
38  2008 JC TIP3P, TIP4Pew, SPC/E LJ Points  Full 
32  2011 KBFF SPC/E LJ Points  Full 
39  2012 MP  LJ Points  Scaled 
40  2014 RDVH SPC/E LJ Points  Full 
41  2014 AH/BK3 BK3 Exponential-6 Gaussian Drude Full 
42  2016 NaCl/ϵ TIP4P/ϵ, SPC/ϵ LJ Points  Scaled 
43  2017 Madrid TIP4P/2005 LJ Points  Scaled 

Activity coefficients measure non-idealities in solution; for electrolytes, these are computed using an infinite-dilution (Henry’s law) reference state. Their determination requires chemical potential (or free energy) values for the electrolyte at a finite concentration, as well as reference-state values. The calculation of free energies or chemical potentials for implicit-solvent ionic systems can be performed directly using Widom test particle insertions44 coupled to molecular dynamics or Monte Carlo sampling of configuration space.45,46 The implicit-solvent approach is not readily transferrable to higher temperatures for which the hydration structure of the solution changes.

Determination of chemical potentials in explicit-solvent simulations of electrolyte solutions is much harder than for implicit-solvent simulations because the strong interactions in such systems result in complex, persistent hydration structures around ions that would be severely disrupted on direct insertion of full-size, full-charge additional entities in solution. Only thermodynamic integration or other such multi-step methods are able to compute the free energies of ions in solution. In addition to the difficulty of obtaining the chemical potential at a finite concentration, the calculation of activity coefficients is impeded by the need to determine the chemical potential at the infinite-dilution reference state, μ. For example, for a 1:1 electrolyte, the formal definition of the mean ionic activity coefficient, γ, at a given concentration m (in moles of solute per kg of solvent) is

2lnγ=βμβμ2lnm,
(1)

where β is the inverse temperature, β = 1/kBT, with kB being Boltzmann’s constant.

Early calculations of ionic activity coefficients in simulations with explicit water include the study of Sanz and Vega47 for a model of NaCl in SPC/E water. Joung and Cheatham10 calculated the mean ionic activity coefficients of aqueous NaCl and KF, as shown in Table III. To avoid the difficulties with estimating the chemical potential of salt at infinite dilution, μ was determined in Ref. 10 by assuming that the mean ionic activity coefficients at the lowest concentration studied, m = 0.038 mol/kg (one ion pair in 1500 water molecules), matched the experimental value.

TABLE III.

Aqueous electrolyte models with known activity coefficients and solubilities.

T = 298.15 KT = 373.15 KT = 473.15 K
ReferenceH2O modelSaltIon modelSsimSexpγsimγexp(m)SsimSexpγsimγexp(m)SsimSexpγsimγexp(m)
10  TIP3P NaCl JC 0.25 1.7 (4)     
  KCl JC 0.21 0.9 (2)     
51  TIP4Pew NaCl JC 0.23 1.8 (6)     
 SPC/E  SD + TF 0.10 3.0 (6)     
 Exp-6  TF 0.01 0.1 (6)     
53  SPC/E NaCl KBFF 0.14 1.0 (6) 0.13 0.5 (6) 0.10 0.4 (5) 
   RDVH 0.93 2.0 (6) 0.78 1.5 (6) 0.64 1.2 (6) 
   JC 0.60 2.3 (6) 0.45 1.6 (6) 0.31 1.1 (6) 
   SD 0.10 1.2 (6) 0.11 0.9 (5) 0.10 0.6 (5) 
  NaI KBFF >1.1 0.6 (12)     
   RDVH 1.1 7.1 (12)     
   JC 0.63 5.8 (12)     
  NaF KBFF 0.01 0.8 (1)     
   RDVH <0.01 0.8 (1)     
   JC 1.3 1.1 (1)     
  KCl KBFF 0.12 1.0 (5)     
   RDVH 0.03 0.7 (5)     
   JC 0.56 1.3 (5)     
65  TIP4P/ϵ NaCl NaCl/ϵ 0.02 0.5 (5) 0.04 0.5 (5) 0.07 0.6 (5) 
 SPC/ϵ  NaCl/ϵ 0.02 0.5 (5) 0.05 0.6 (5) 0.10 0.7 (5) 
66  SWM4 NaCl AH/SWM4 0.13 0.5 (3)     
 BK3  AH/BK3 0.16 1.0 (6)     
52  BK3 NaCl AH/BK3 0.10 1.0 (4) 0.11 0.8 (4) 0.10 0.7 (4) 
43  TIP4P/2005 NaCl Madrid 0.93 1.3 (6)     
   JC 0.57 1.4 (6)     
T = 298.15 KT = 373.15 KT = 473.15 K
ReferenceH2O modelSaltIon modelSsimSexpγsimγexp(m)SsimSexpγsimγexp(m)SsimSexpγsimγexp(m)
10  TIP3P NaCl JC 0.25 1.7 (4)     
  KCl JC 0.21 0.9 (2)     
51  TIP4Pew NaCl JC 0.23 1.8 (6)     
 SPC/E  SD + TF 0.10 3.0 (6)     
 Exp-6  TF 0.01 0.1 (6)     
53  SPC/E NaCl KBFF 0.14 1.0 (6) 0.13 0.5 (6) 0.10 0.4 (5) 
   RDVH 0.93 2.0 (6) 0.78 1.5 (6) 0.64 1.2 (6) 
   JC 0.60 2.3 (6) 0.45 1.6 (6) 0.31 1.1 (6) 
   SD 0.10 1.2 (6) 0.11 0.9 (5) 0.10 0.6 (5) 
  NaI KBFF >1.1 0.6 (12)     
   RDVH 1.1 7.1 (12)     
   JC 0.63 5.8 (12)     
  NaF KBFF 0.01 0.8 (1)     
   RDVH <0.01 0.8 (1)     
   JC 1.3 1.1 (1)     
  KCl KBFF 0.12 1.0 (5)     
   RDVH 0.03 0.7 (5)     
   JC 0.56 1.3 (5)     
65  TIP4P/ϵ NaCl NaCl/ϵ 0.02 0.5 (5) 0.04 0.5 (5) 0.07 0.6 (5) 
 SPC/ϵ  NaCl/ϵ 0.02 0.5 (5) 0.05 0.6 (5) 0.10 0.7 (5) 
66  SWM4 NaCl AH/SWM4 0.13 0.5 (3)     
 BK3  AH/BK3 0.16 1.0 (6)     
52  BK3 NaCl AH/BK3 0.10 1.0 (4) 0.11 0.8 (4) 0.10 0.7 (4) 
43  TIP4P/2005 NaCl Madrid 0.93 1.3 (6)     
   JC 0.57 1.4 (6)     

Efficient methods for obtaining chemical potentials and activity coefficients of electrolyte solutions from simulations have been developed in recent years. These methods also allow the determination of solubilities, when combined with methods to obtain the chemical potentials of the salt crystals.48 Moučka et al.9,49,50 calculated chemical potentials (from which activity coefficient expressions were also obtained) and solubilities for aqueous NaCl using the SPC/E water model24 with several different ion models. They used the “Osmotic Ensemble Monte Carlo” (OEMC) method of Moučka et al.,49 which is an open-ensemble method involving gradual insertion of ions to maintain an input value of the chemical potential.

An approach developed in the author’s group51,52 is based on optimized gradual insertions of neutral ion groups (pairs for 1:1 salts) into a solution of given concentration to determine the chemical potential in the fluid phase. The total sampling time required to obtain reasonable statistics for the calculation of activity coefficients is on the order of μs, a time scale that currently can only be accessed for computationally efficient classical force fields. Extrapolation using Debye–Hückel theory to the infinite-dilution limit is used to estimate the reference chemical potential needed for the activity coefficient calculation. The methodology has been further refined53 by utilizing the Davies equation54 for the extrapolation of the excess chemical potential to infinite dilution. The main difference of this approach relative to the OEMC approach is the flipping of the independent variable from the chemical potential to its conjugate, concentration. The results for the chemical potentials as a function of concentration from the two methods have been confirmed to be in agreement within their respective uncertainties.12 

One alternative to thermodynamic integration methods for obtaining the chemical potential of the salt is to compute the chemical potential of the solvent and then use the Gibbs–Duhem relationships between the chemical potentials of salt and solvent to obtain the latter, as first proposed in Ref. 50. In addition to standard methods involving insertions, this can also be done by direct simulations of the osmotic pressure.55 This “osmotic pressure for activity of solvents” (OPAS) method relies on a virtual semipermeable membrane to confine the electrolyte solution and performs the measurements of the force acting on the membrane. The osmotic pressure is an important property in itself.56 Alternatively, the density of the vapor phase (of pure solvent) at equilibrium with the liquid solution can be obtained and used to compute the activity coefficient of the electrolyte from the Gibbs–Duhem relationships.57 This approach has been used to obtain activity coefficients of aqueous dysprosium nitrate solutions.58 

A fundamental challenge in the simulations of electrolyte solutions is that Coulombic interactions are of infinite range, decaying only as the inverse power of the distance, r−1. At first glance, this would imply that, in a uniform system, “far away” interactions are more important than “nearby” ones, given that the number of interaction centers at distance r grows as 4πr2. This is the reason for which systems dominated by gravitational interactions (which also decay as r−1) do not have intensive properties and cannot be described by equilibrium thermodynamics. Of course, ionic systems, unlike gravitational ones, do have thermodynamic properties, because of electroneutrality: at long distances, there is charge screening and effective interactions decay exponentially.

Treatment of the long-range Coulombic interactions is a significant issue in simulations of electrolyte solutions. Modern simulations of such systems are almost always performed using Ewald sums59,60 or equivalent methods, rather than crude truncation methods used in early simulations. In the Ewald approach, the summation of the electrostatic interactions is performed by taking advantage of the periodic character of the simulated system. In the simulations of solutions, there is typically a non-zero dipole within the unit cell, which, in turn, interacts with its copies in the periodically replicated boxes and with an external dielectric medium surrounding the macroscopic collection of replicas at infinite distance. The medium at infinity can be vacuum (ϵ = 1), a conductor (ϵ = ∞), or a medium with an intermediate relative permittivity.

There are well-documented finite-size effects on the chemical potential of ions obtained through simulations.61 Of direct relevance to the calculation of chemical potentials and activity coefficients in electrolyte solutions is the study of Thompson and Sanchez,62 who obtained the system-size dependence of free energy for “weakly coupled” ionic fluids using theoretical arguments, and the recent simulation study from the author’s group63 who used implicit-solvent simulations to probe the magnitude of finite-size effects on the free energies. Figure 1 illustrates the magnitude of these effects at a low (top panel) and intermediate (bottom panel) concentration of a model 1:1 electrolyte. At low concentrations, because there are only a few ion pairs in the simulation box, these effects are very strong, while ion screening at higher concentrations reduces their magnitude. The different curves correspond to different Ewald boundary conditions at infinite distance and different simulation methods, specifically thermodynamic integration coupled with molecular dynamics, and direct insertions of ion pairs in grand canonical Monte Carlo. The curves can all be extrapolated to a common value corresponding to the “true” chemical potential for the physical system of interest. These finite-size effects, if left uncorrected, lead to systematic errors of over 0.1 units in the logarithm of the activity coefficients (or 10% in the activity coefficient values).63 Clearly, a careful examination of systems’ size effects and extrapolation to infinite system size are necessary for obtaining accurate results for free energies and activity coefficients in simulations of electrolyte solutions.

FIG. 1.

Dependence of the chemical potential on system size for a model of NaCl solutions at concentrations of 0.0066 mol/l (top panel) and 0.14 mol/l (bottom panel) from Ref. 63. BC stands for boundary conditions at infinite distance in the Ewald summations. Error bars and shading indicate statistical uncertainties for the individual simulations and the fitted (dashed) lines, respectively. Reprinted with permission from J. M. Young and A. Z. Panagiotopoulos, J. Phys. Chem. B 122, 3330 (2018). Copyright 2018 American Chemical Society.

FIG. 1.

Dependence of the chemical potential on system size for a model of NaCl solutions at concentrations of 0.0066 mol/l (top panel) and 0.14 mol/l (bottom panel) from Ref. 63. BC stands for boundary conditions at infinite distance in the Ewald summations. Error bars and shading indicate statistical uncertainties for the individual simulations and the fitted (dashed) lines, respectively. Reprinted with permission from J. M. Young and A. Z. Panagiotopoulos, J. Phys. Chem. B 122, 3330 (2018). Copyright 2018 American Chemical Society.

Close modal

Typical results for the activity coefficients of NaCl using different force fields are shown in Fig. 2, taken from Ref. 52. Experimental data are shown as lines and simulation results as points. According to the Debye–Hückel limiting law, the logarithm of the mean ionic activity coefficient, ln γ, follows a straight line at low concentrations with a slope that depends only on the dielectric permittivity of the solvent and the charge of the ions, when plotted against m. The physical reason for the decrease in the activity coefficients at low concentrations is the mutual attraction of the oppositely charged ions. Simulation results plotted in Fig. 2 closely follow the Debye–Hückel limiting law at low concentrations, since the relative permittivity for the chosen water models is close to the experimental value and full charges are assumed for the ion models selected. At higher concentrations, the activity coefficients flatten off and then increase in magnitude, corresponding to less favorable ion–ion interactions because of the changing structure of the electrolyte solutions, the presence of repulsive ion cores, and crowding in the hydration shells of individual ions. At room temperature (red-colored symbols and line), the BK3 polarizable model, shown as filled symbols, does remarkably well over the whole concentration range studied, especially taking into account that the activity coefficients were not used in the development of the model parameters for BK3. The SPC/E + SD water and ion model combination (open symbols) does less well, systematically overestimating the activity coefficients at high concentrations. At the highest temperature shown (473.15 K, black-colored symbols and line), both polarizable and non-polarizable models significantly underestimate the activity coefficients, resulting in more favorable ion–ion interactions relative to the experimental observations.

FIG. 2.

Activity coefficients of NaCl for the selected ion/water models as a function of temperature from Ref. 52. Reprinted with permission from Jiang et al., J. Chem. Theory Comput. 11, 3802 (2015). Copyright 2015 American Chemical Society.

FIG. 2.

Activity coefficients of NaCl for the selected ion/water models as a function of temperature from Ref. 52. Reprinted with permission from Jiang et al., J. Chem. Theory Comput. 11, 3802 (2015). Copyright 2015 American Chemical Society.

Close modal

A broader set of results for the activity coefficients are summarized in Table III for models for which the solubility values have also been reported. This table lists the ratio of the simulation to the experimental values for the activity coefficients, γsim/γexp, at a specific concentration m, in units of mol/kg of solvent. Deviations of calculated activity coefficients from experimental results generally increase at higher concentrations (as shown in Fig. 2), so the activity coefficient ratios in this table are reported at the highest concentration studied in the corresponding reference. Citations for the water and ion models are listed in Tables I and II, respectively. As shown in Table III, large errors in the activity coefficients are quite common for the majority of the model combinations and conditions studied.

A large number of additional activity coefficients for alkali halides are available from Ref. 64 for the SPC/E water24 plus RDVH40 ion models, but since the corresponding solubilities were not reported, these results are not shown in Table III. Reference 9 has presented results for the solubilities and chemical potentials of 13 aqueous NaCl force fields in SPC/E water, including the SD and JC force fields later studied in Ref. 53. Even though it may be possible to reanalyze the reported chemical potential data into activity coefficients, this would require assumptions about the reference state chemical potential, so this analysis was not attempted here.

Some trends are apparent in this table with respect to the activity coefficients:

  • Most models overestimate the experimental values of activity coefficients at room temperature, except for the KBFF models (which were fitted to these values) and the BK3 models (which were not). It should be noted here that for many of the models studied, the activity coefficient comparison in Table III is done at concentrations above the solubility limit for the corresponding model; this is not a significant issue for most cases, as long as homogeneous nucleation does not occur over the time scale for the measurement of the activity coefficients.

  • Activity coefficients for practically all model combinations studied decrease at a faster rate with temperature than the rate of change seen experimentally; as a result, models that do well at room temperature underestimate the activity coefficients at high temperatures (e.g., BK3), and some models that were not accurate at room temperatures become better at some elevated temperature (e.g., SPC/E + SD at 373 K, or SPC/E + JC or RDVH at 473 K). This, however, is little more than fortuitous cancellation of errors. Unfortunately, all the available data for the temperature dependence of the activity coefficients are for a single salt, namely, NaCl. Thus, it is not clear if the trends shown in Fig. 2 are generally valid for other salts.

  • The Madrid scaled-charge model for NaCl does relatively well for the activity coefficients at room temperature. However, the model overpredicts the experimental γ values over the whole concentration range. This stands in contrast to the full-charge models, which generally fit the activity coefficients at low-to-intermediate concentrations quite well. Activity coefficients at elevated temperatures have not been determined for the Madrid scaled-charge model.

There are relatively few results for activity coefficients of multivalent salts in explicit water. Mamatkulov et al.67 calculated activity coefficient derivatives for alkaline earth chlorides in order to parameterize force fields for the divalent cations for use with the SPC/E water model. Not surprisingly, the resulting force fields do well for the activity coefficients of CaCl2 at room temperature, but start deviating from the experimental results at higher temperatures,68 in the same direction of underprediction of the activity coefficients as the 1:1 electrolyte results shown in Fig. 2. The force field of Deublein et al.69 for CaCl2 also does well at room temperature, but overpredicts significantly the activity coefficients at high temperatures.68 Finally, a reduced-charge model70 does poorly for the activity coefficients because it does not follow the Debye–Hückel limiting law at low concentrations.

There are two main approaches used to determine solubilities by simulation, namely, the indirect free energy calculation route and the direct coexistence method. To obtain the solubility using the free energy route, the chemical potential of the solid phase is determined, typically using the Einstein crystal method.48 The dependence of the chemical potential on concentration is calculated in the solution phase, often from the free energy change of an ion pair insertion, using methods similar to those described in Sec. III A. The concentration at which the chemical potential of the solution is equal to the chemical potential of the solid crystal is the solubility. In this context, it is of interest to note that Li et al.71 presented “a systematic extension of the Einstein crystal method to calculate the absolute solid free energies of molecular crystals at arbitrary temperatures and pressures” and combined it with a cavity method to compute free energies in solution in order to obtain the solubilities of sparingly soluble solutes. This method was applied, among other systems, to the sparingly soluble CaCO3 salts calcite and aragonite.72 In other work involving chemical reactions, Noroozi and Smith modeled equilibrium constants and the reactive absorption of CO2 in aqueous amine solutions involving multiple electrolytes and reactions.73,74 The indirect free energy calculation approach requires multiple simulations, frequently at “uninteresting” state points, which are computationally expensive and complex. However, the solution chemical potential is also required for the activity coefficients, and the resulting solubility values are more reliable than the solubility from the direct coexistence method for salts.75 

In the direct coexistence method, a solid slab is placed next to a solution slab. The system is allowed to reach a steady state, which requires μs-long molecular dynamics calculations for salt systems77 because of the slow rates of ionic crystal dissolution or formation. The concentration in the bulk of the solution phase after equilibration is taken as the solubility. The direct coexistence method requires no free energy calculations but suffers from strong system-size effects.76 The reason for these effects is the preferential adsorption of ions from solution to the crystal interfaces that lead to surface dipoles interacting through the low-dielectric crystal slabs.76 This, in turn, leads to the strong dependence of the solubility values on slab thickness. This is illustrated in Fig. 3, reproduced from Ref. 76, which shows the apparent solubility of the JC model for NaCl38 in SPC/E water24 at 298 K, using both direct coexistence and free energy methods. Direct coexistence simulations with different crystal slab aspect ratios are shown as different colored symbols as a function of the width of the crystal slab or crystallite. The values computed from chemical potential calculations are shown as a horizontal line with associated uncertainty, since in this method, the simulated system is the bulk crystal under periodic boundary conditions in all directions. There is a strong system-size dependence for the direct coexistence calculations. Early simulations were performed with relatively narrow slabs (right-hand side of Fig. 3), leading many researchers to conclude that the solubility for this model combination is near 6 mol/kg, close to the experimental value. However, for sufficiently large slabs, the results from the direct coexistence calculations are significantly lower, agreeing within simulation uncertainties with the value 3.7 ± 0.2 determined from free energy calculations.

FIG. 3.

Apparent solubility of the JC model for NaCl in SPC/E water at 298 K as a function of the inverse width of the crystal slab for direct coexistence calculations (red, green, and purple symbols—different letters indicate different aspect ratios of the simulated systems), a spherical crystallite (SC, blue circle), and from free energy calculations (black horizontal line with uncertainties marked by dashed horizontal lines). Reproduced with permission from Espinosa et al., J. Chem. Phys. 145, 154111 (2016). Copyright 2016 AIP Publishing LLC.

FIG. 3.

Apparent solubility of the JC model for NaCl in SPC/E water at 298 K as a function of the inverse width of the crystal slab for direct coexistence calculations (red, green, and purple symbols—different letters indicate different aspect ratios of the simulated systems), a spherical crystallite (SC, blue circle), and from free energy calculations (black horizontal line with uncertainties marked by dashed horizontal lines). Reproduced with permission from Espinosa et al., J. Chem. Phys. 145, 154111 (2016). Copyright 2016 AIP Publishing LLC.

Close modal

In addition to the activity coefficients mentioned earlier, Table III collects information on solubilities of aqueous electrolyte models for which both properties are directly available from simulations. This table lists the system and reference in which the results were reported, the corresponding water and electrolyte models (from Tables I and II), and the ratio of the simulation to the experimental values of the solubility, Ssim/Sexp. Additional solubilities have been reported in the literature, e.g., for the polarizable MAH/BK3 model,78 which is a modification of the AH/BK3 model that gives a much improved solubility prediction at room temperature. For non-polarizable models, a recent study79 provides a comprehensive summary of solubility values reported in the literature and provides new direct-coexistence results for many existing models as well as a re-optimized set of LJ model parameters for NaCl and KCl that reproduce the experimental solubility at room temperature.

Several trends for the solubilities are worth pointing out:

  • All model combinations for NaCl underpredict the experimental solubilities, sometimes severely so, at room temperature as well as elevated temperatures. In this respect, some models are particularly poor, e.g., the Exp-6 model of Errington and Panagiotopoulos,25 which was optimized to the phase behavior of pure water, gives an equilibrium solubility of only 1% of the experimental value when combined with the TF NaCl model with which it shares a functional form for the repulsive interactions. Similarly, the NaCl/ϵ model42 gives a solubility of 2% of the experimental value at room temperature, with the results improving somewhat at higher temperatures. Such models with very low equilibrium solubilities are likely to lead to rapid aggregation and precipitation for simulations of salt solutions near the experimental solubility limit, leading to invalid results. The underprediction of solubilities also extends to the polarizable models studied. The fact that solubility values common electrolyte/water model combinations are too low has also been pointed out in many prior articles, e.g., in Ref. 13.

  • The situation for other alkali halides is a bit more mixed, with preponderance of underestimated solubilities, but some model combinations getting close (or even above) the experimental values. Once more, there are model combinations with equilibrium solubilities only a few percent of the experimental values.

  • The temperature dependence of the solubility is also not correctly predicted by most models, with the ratio of Ssim/Sexp generally becoming smaller as temperature is increased. Since the solubilities at room temperature are typically underestimated, this means that deviations from experiments become greater at higher temperatures.

  • The Madrid scaled-charge model for NaCl does well for the solubility at room temperature, as it also did for the activity coefficients at high concentrations. The temperature dependence of the solubility is not currently known for this model. A number of additional scaled-charged models are available for different alkali and alkaline earth halides,80 but their solubilities and activity coefficients have yet to be determined.

There is a strong inverse correlation between the ability of a model to predict correctly experimental solubilities and activity coefficients. Specifically, model combinations that do well for the activity coefficients (e.g., TIP3P + JC for KCl at 298 K) do poorly for the solubilities, and the converse is also true. The fact is that the deviations from experimental behavior are systematic, and the correlations between solubility and activity coefficient predictions strongly point to “missing physics” in the models. These may be the result of fundamental limitations of the existing modeling framework that may not be addressable through simple tweaks in the parameters. For the fixed-point-charge models of Table III, it could be argued that the “missing physics” is the polarizability of the ions and water. However, polarizable models show similar trends, perhaps to a somewhat lesser extent—this deficiency could, however, be the consequence of the inadequate optimization of the polarizable model parameters. This last argument is harder to make for the non-polarizable models for which parameter optimization has been going on for a long time. This hypothesis is bolstered by the very recent study of Dočkal et al.,81 who obtained a new parameterization for a NaCl polarizable force field starting from the AH/BK3 and MAH/BK3 models, which does remarkably well for both crystal and solution chemical potentials, leading to solubility predictions near the experimental values.

For multivalent salts, the only system for which some information is known for both activity coefficients and solubilities is CaCl2. Moučka et al.82 were the first to report the solubility using the Mamatkulov et al. force field67 in combination with SPC/E water. The room-temperature hexahydrate phase was found to be unstable, while the solubility of the high-temperature dihydrate phase was found to be lower than the experimental value by around 50%. Based on these results, the authors concluded that the force field “does not provide a good description of the experimental properties considered.” The same model was studied by the author’s group,68 along with another full-charge model69 and a reduced-charge model.70 All models were found to underpredict the experimental solubilities of the various hydrate phases of CaCl2. However, extremely slow dynamics at high electrolyte concentrations made equilibration of the systems problematic, and the reported values may not correspond to thermodynamic equilibrium. Because of these sampling limitations, these results are not given in Table III.

The transport properties of electrolyte solutions are of interest for applications and for the fundamental understanding of the dynamic behavior of such systems. These properties include diffusivities of the ions and water, solution conductivities, and viscosities. High-quality experimental data for these properties are available for many single salts in water83,84 as well as for some mixtures, usually at room temperature. Transport coefficients can be obtained from equilibrium molecular dynamics simulations using the Green–Kubo transport relationships or from non-equilibrium simulations in which gradients are imposed. The former (equilibrium) approach is preferred because it does not require extrapolations to zero gradient conditions and is quite feasible using modern codes and computational facilities.

Ionic mobilities are directly related to solution conductivity and are thus important for electrochemical applications. Koneshan and co-workers85,86 performed early simulation studies of ion mobilities in the bulk, suggesting the presence of distinct maxima with respect to ionic radii. Joung and Cheatham10 obtained ion diffusion coefficients for a number of water and ion model combinations.

In addition to their intrinsic value, the dependence of the transport coefficients on concentration conveys significant information on the structure of the solution, as modified by the presence of ions. For example, the viscosity η of electrolyte solutions at relatively low concentrations (where ion–ion interactions are not important) is usually well-represented by the following expression:84 

η/η0=1+Am+Bm.
(2)

In this expression, η0 is the viscosity at zero concentration. Coefficient A can be derived theoretically from the drag that the ion atmosphere imposes on the solution; its value depends on the viscosity of the solvent, its permittivity, and the ion mobilities and valences. The “Jones–Dole” coefficient B in the viscosity expression is connected to the “degree of hydration” of the ions.87 The sign of the B coefficient is connected to the Hofmeister series characterization of ions as “cosmotropes” (structure-makers, positive B coefficients) or “chaotropes” (structure-breakers, negative B coefficients). Until very recently, there had been no simulation data for electrolyte solution viscosities of sufficiently high quality to allow for direct comparisons of the A and B coefficients for specific molecular models to experimental data. For example, Orozco et al.88 obtained viscosities at intermediate and high concentrations for a number of non-polarizable water + ion models, and Jiang et al.52 obtained the shear viscosity of aqueous NaCl for the BK3 + AH/BK3 polarizable model. Neither study attempted to extract the viscosity coefficients A and B from the results in order to compare them to theoretical and experimental values.

A landmark study of the diffusivity and viscosity of aqueous electrolyte solutions89 revealed that none of the common fixed-point-charge models predict the correct qualitative effects on solvent mobility. Contrary to experimental results, both “structure making” and “structure breaking” ions were observed to decrease the mobility of water in the simulations, resulting in lower diffusivities and higher viscosities. Even the two polarizable models studied in Ref. 89 did not perform well against experimental data. The scaled-charge multi-body potential of Kann and Skinner34 does better but still underestimates the magnitude of the enhanced diffusivity of water for CsI, a structure-breaking salt.

Recently, Yue and Panagiotopoulos90 performed a systematic study of the viscosity B-coefficients and water mobility for several of the ion models listed in Table II in combination with the SPC/E, TIP4P/2005, and BK3 water models of Table I. The results from this study for the viscosity B-coefficients of NaCl solutions are presented in Fig. 4. Several other alkali halide salts (KCl, NaI, NaBr, LiCl, and NaF) were studied in Ref. 90, with comparable results. In order to obtain the transport properties of sufficiently high quality to allow for the determination of the B-coefficients, 100 ns long molecular dynamics simulations were performed for 100 independent initial configurations so that the total sampling time was 10 µs at each of six concentrations for each model studied.

FIG. 4.

The linearized form of the viscosity Jones–Dole equation for NaCl solutions at various concentrations from Ref. 90. The slope of the lines corresponds to the B-coefficient of Eq. (2), with averages (dashed lines) and statistical uncertainties (shaded regions) given as 95% confidence intervals. Models used are SPC/E + JC (“nonpolarizable”), TIP4P/2005 + MP (“scaled-charge”), and BK3 + AH/BK3 (“polarizable”), as listed in Tables I and II. Experimental values (solid line) are from Ref. 84. Reprinted with permission from S. Yue and A. Z. Panagiotopoulos, Mol. Phys. 117, 3538 (2019). Copyright 2019 Taylor and Francis.

FIG. 4.

The linearized form of the viscosity Jones–Dole equation for NaCl solutions at various concentrations from Ref. 90. The slope of the lines corresponds to the B-coefficient of Eq. (2), with averages (dashed lines) and statistical uncertainties (shaded regions) given as 95% confidence intervals. Models used are SPC/E + JC (“nonpolarizable”), TIP4P/2005 + MP (“scaled-charge”), and BK3 + AH/BK3 (“polarizable”), as listed in Tables I and II. Experimental values (solid line) are from Ref. 84. Reprinted with permission from S. Yue and A. Z. Panagiotopoulos, Mol. Phys. 117, 3538 (2019). Copyright 2019 Taylor and Francis.

Close modal

Results from Refs. 89 and 90 clearly illustrate that non-polarizable, fixed-charge ion and water models do not represent accurately the transport properties of aqueous electrolyte solutions at low concentrations, with all ions, regardless of their Hofmeister classification, showing kosmotropic behavior. The viscosity B-coefficients obtained are higher than the experimental values, even though the relative order of different anions and cations is qualitatively correct.90 Scaled charge models are significantly better in this respect, with the predicted dependence of viscosity on concentration being closer to experiments. The BK3 + AH/BK3 model combination studied in Ref. 90 and the newly proposed Dočkal et al.81 polarizable force field are similar in performance for the transport coefficient to those of scaled-charge models.

Recent studies of a coarse-grained solution model of Andreev et al.91,92 are relevant in the context of transport properties for electrolyte solutions. Their model consists of charged LJ spheres (ions) in an LJ solvent representing water. By setting the strength of the ion–solvent attractive interactions using the measured enthalpies of solvation of ions in water, Andreev et al. were able to reproduce trends of diffusivity for both structure-making and structure-breaking salts. They also examined the connection between the B coefficient and the partial molar volumes of the salts. This “brutally coarse-grained model of electrolyte solutions” thus captures effects of salts on solution transport properties that are missed by more realistic atomistic models, but at the cost of being unable to describe higher electrolyte concentrations, salt activities and solubilities, or nucleation rates, which are the subject of Sec. VI.

The rates and mechanisms of precipitation from solution to solid (crystal) phases are of fundamental scientific interest, with important ramifications in many disciplines, e.g., catalysis, petroleum and pharmaceutical industries, geochemistry, atmospheric sciences, and biology. There is, however, a significant gap in our understanding of electrolyte precipitation as the highly non-equilibrium nature of the nucleation process makes probing the microscopic mechanism of nucleation extremely difficult.93 

Nucleation is a rare event, except at extremely high supersaturations, where spontaneous nucleation may be observed within a few hundred nanoseconds, time scales accessible in standard molecular dynamics calculations.94,95 Determination of nucleation rates by simulation at more realistic supersaturations relies either on free energy calculations or on a variety of path-sampling methods.96 Free energy calculations typically use umbrella sampling97 or metadynamics98 to allow a system to move over the nucleation barrier. The results from such calculations can be quite sensitive to the choice of order parameter used to describe the transition, and the calculation of actual nucleation rates requires additional (kinetic) information.

An alternative to free energy calculations is provided by path-sampling methods,99 which do not perturb the natural dynamics of a system and are independent of assumptions on the order parameter relevant for the underlying transition. A particularly promising path-sampling approach is forward-flux sampling.100,101 The method relies upon launching a large ensemble of molecular dynamics simulations from specific “milestones” along the nucleation process and counting the fraction of members of these ensembles that reach the next milestone. The milestones can be defined through, e.g., the size of the largest crystal present in solution, but the resulting nucleation rates do not depend on this choice. The approach has been comprehensively reviewed in this journal in a recent perspective article by Hussain and Haji-Akbari.102 

Path-sampling simulations are computationally very demanding, as they require following a large number of independent trajectories for sufficiently long times so that significant structural rearrangements can occur, in thermodynamic parameter regions where dynamics can be quite slow. Thus, it is appealing to look for faster alternatives. Seeding simulations77 use pre-formed crystals, usually of spherical shape, embedded in a metastable liquid. The growth or dissolution of the seed is observed, and the critical nucleus size determined when the crystal is stable with time. Nucleation rates are then extracted using the equations of classical nucleation theory. The nucleation rate is a strong function of supersaturation, the ratio of solute concentration to the equilibrium solubility. Thus, it is of critical importance to have a good estimate of the equilibrium solubility for the water and ion model combination being studied prior to attempting a calculation of nucleation rates. For example, the initial study of Zimmermann et al.103 used an earlier direct-coexistence estimate of the solubility of JC NaCl in SPC/E water, leading to errors of many decades in order of magnitude for the nucleation rate estimate.104 Even when the correct solubility is used, it turns out that for salt crystallization, nucleation rates calculated from the seeding method depend strongly on the order parameter used to identify the extent of the crystalline region in the simulation box.104 While seeding simulations are significantly less computationally demanding than forward-flux sampling methods, it seems that reliable results can only be obtained by calibrating the selection of optimal order parameters and clustering rules against more rigorous methods.

Jiang et al.105 used forward-flux sampling calculations to obtain the first rigorous estimate of the homogeneous nucleation rate for the JC NaCl force field combined with the SPC/E water model at 298 K and 1 bar. The results, shown as blue up-triangles in Fig. 5, suggest that this model combination significantly underestimates the nucleation rate at all relative supersaturations for which experimental data106 (red crosses) are available. Note the highly compressed logarithmic ordinate scale of this figure. A subsequent study107 extended the calculations to the polarizable MAH/BK3 model78 and the Madrid reduced-charge model,43 which reproduces the experimental equilibrium solubility quite well, as indicated in Table III. The reduced-charge model completely fails with respect to the nucleation rate, resulting in underprediction of the experimentally observed rates by 20 orders of magnitude. This is perhaps not surprising, given that the properties of the crystal phase and the solution–crystal interface are not expected to be described well by the reduced-charge model. Interestingly, the polarizable MAH/BK3 model shows a qualitatively different dependence of the nucleation rate on supersaturation relative to the two non-polarizable models studied. Figure 3 of Ref. 107 (not shown here) indicates that the underlying cause for this difference is that the melt-solution interfacial tension is lower and closer to the experimental value for the polarizable model. Even though this is only a single study of three specific models for NaCl, the results are intriguing: if they can be shown to be valid for other systems, the implications are that polarizable models make a significant difference with respect to the description of interfacial free energies and nucleation rates, despite their limited success thus far with respect to solubilities and activity coefficients.

FIG. 5.

Homogeneous nucleation rate J for aqueous NaCl as a function of the supersaturation ratio S, from Ref. 105. Blue up-triangles are for the JC NaCl + SPC/E water model combination, red crosses are experimental data from Ref. 106, squares are for the MAH/BK3 polarizable model combination of Ref. 78, and down-triangles are for the Madrid reduced-charge model. Reprinted with permission from Jiang et al., J. Chem. Phys. 148, 044505 (2018). Copyright 2018 AIP Publishing LLC.

FIG. 5.

Homogeneous nucleation rate J for aqueous NaCl as a function of the supersaturation ratio S, from Ref. 105. Blue up-triangles are for the JC NaCl + SPC/E water model combination, red crosses are experimental data from Ref. 106, squares are for the MAH/BK3 polarizable model combination of Ref. 78, and down-triangles are for the Madrid reduced-charge model. Reprinted with permission from Jiang et al., J. Chem. Phys. 148, 044505 (2018). Copyright 2018 AIP Publishing LLC.

Close modal

The mechanism of nucleation in salt solutions is an important problem on which significant progress has been made in recent years. Chakraborty and co-workers94,95,108 studied nucleation of NaCl crystals from highly supersaturated solutions and observed a two-step mechanism108 in which an ion “patch” initially forms within which a second spatial ordering step takes place to produce the salt crystals. This concept was confirmed recently by Jiang et al., 109 who also determined the chemical potential values of the supersaturated solution as a function of concentration. As shown in Ref. 109, the two-step mechanism, first to an “amorphous salt” phase and then to the crystal takes place at the concentration that corresponds to the spinodal (limit of thermodynamic stability) for the solution, a phenomenon, which has also been observed earlier.9,12,49 At lower concentrations, a simple, one-step nucleation-and-growth mechanism dominates. A similar initial transition to an amorphous phase upon crossing the spinodal has been observed in crystallization from a LJ mixture.110 Progress is also being made on elucidating the ion dissolution mechanism and kinetics111,112 and effects of salts on the nucleation rate of ice.113 

This article provides an overview of simulation methods, molecular models, and results related to the collective properties of electrolyte solutions, specifically activity coefficients, solubilities, transport properties, and nucleation rates. Algorithmic developments in recent years, coupled with the increased availability of high-performance computing facilities, have enabled calculations that were impractical a decade or so ago. Even now, however, sophisticated models properly describing multi-body interactions in electrolyte solutions pose significant computational challenges for properties that require large systems and long sampling times. Fully ab initio methods are still many orders-of-magnitude too slow for practical computations of many of the collective properties discussed in the present review. For this reason, the models discussed here are relatively simple, representing water and ions using a few interaction centers. Many of the models only include fixed point charges for the electrostatic interactions.

The vast majority of molecular models for water and ions in solution have not been parameterized for the properties we have examined here. Nevertheless, many of them do reasonably well with respect to modeling of many aspects of the collective behavior of electrolyte solutions. This is a consequence of utilizing physically realistic models that capture the main interactions in a system and is a key advantage of molecular simulation over phenomenological equations that can only describe a narrow set of properties over the range of conditions for which experimental data have been fitted to them. However, when examining the quantitative performance and transferability of the available models, significant limitations emerge. Specifically, many models underpredict the solubility of salts in water at room temperature. There is an inverse correlation between the quality of representation of the activity coefficients over the complete composition range and that of the solubility predictions. Moreover, the temperature dependence of both solubility and activity coefficients is incorrect. The transport properties of the solutions are described more accurately by reduced-charge and polarizable models relative to non-polarizable, full-charge models.

A notable gap in the literature exists for the collective properties of interest to the present article with respect to the range of salt types for which data are available. By far, most studies of activity coefficients and solubilities are for NaCl; few studies examine the temperature dependence of these properties. It would be of significant interest for future work to examine how the trends identified for NaCl carry over to different salts, including ion combinations with different classifications with respect to chaotropic/kosmotropic behavior.

Charge scaling is an excellent way to incorporate approximately some of the effects of polarizability that real ions and water possess, at a much lower computational cost than full polarizability. However, as stated by Zeron et al.,80 “One should recognize from the very beginning that the scaling prevents an accurate description of the solid phase and/or the molten salts. The lack of transferability to other phases is the price to pay to improve the description of the solutions.” Force fields based on scaled charges can provide good thermodynamic properties and have clear advantages over full-charge, non-polarizable models for the transport properties of solutions, but they cannot represent accurately the activity coefficients at low and intermediate concentrations. Such models also do not perform well for the properties of the crystal phases and for the prediction of nucleation rates. A final issue arises when chemical transformations among charged species in molten salts are considered in systems relevant for fuel cells and carbon sequestration applications. Because of the charge scaling, the free energies of transfer from the gas phase are incorrect and it is not possible to perform accurate calculations of chemical equilibria.

Fully polarizable, many-body force fields are starting to emerge that appear to be promising for the description of the properties of interest to this article. Nevertheless, most existing polarizable models are not able to represent accurately both activity coefficients and solubilities and do not have the correct temperature dependence for these properties, with the notable exception being the Dočkal et al. force field.81 At least to some extent, these deficiencies may be the result of inadequate parameterization of the models, which have not “been around” for as long as the simpler non-polarizable versions. For many recent polarizable models, e.g., the E3B (explicit 3-body) model of Tainter and co-workers114 or the HBP model for water of Jiang et al.,115 compatible, verified ion models have yet to be developed. Polarizable models seem to have a distinct advantage over non-polarizable ones for the description of nucleation rates of salts from solution, but this conclusion is currently based on a limited set of results for a few models and needs to be verified with respect to its broad applicability.

It seems clear that in order to describe properly charge-transfer and other effects that seem to be important for a quantitative description of the collective properties of aqueous electrolyte solutions, we will eventually need to take into account explicitly their quantum chemical character. Currently, ab initio calculations are still many orders of magnitude too slow—for example, first-principles molecular dynamics methods,116 which have the same accuracy as density functional theory, are currently limited to system sizes of only hundreds of atoms and time scales of less than 100 ps.117 One way to extend the time and length scales that can be handled with quantum chemical methods is to use machine learning to obtain representations of the potential energy surface of the underlying system.118,119 The “deep potential” approach of Han and co-workers117 seems particularly promising in this respect. It would be of significant interest to determine if modern functionals such as SCAN120 could be used to describe the collective properties of electrolyte solutions at a fundamentally improved level relative to empirical force fields.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Financial support for work on electrolyte solutions in the author’s group was provided by the Office of Basic Energy Sciences, U.S. Department of Energy, most recently under Award No. DE-SC0002128. Significant computational resources were provided by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High-Performance Computing Center and Visualization Laboratory at Princeton University. I would like to thank Professor Baron Peters for providing helpful references and the anonymous reviewers for the careful reading of this manuscript and their excellent suggestions.

1.
P. W.
Atkins
,
Physical Chemistry
, 4th ed. (
W. H. Freeman and Co.
,
NY
,
1990
).
2.
J. W.
Tester
and
M.
Modell
,
Thermodynamics and Its Applications
, 3rd ed. (
Prentice Hall
,
Upper Saddle River, NJ
,
1996
).
3.
J. C.
Rasaiah
,
D. N.
Card
, and
J. P.
Valleau
,
J. Chem. Phys.
56
,
248
(
1972
).
4.
S. J.
Marrink
,
H. J.
Risselada
,
S.
Yefimov
,
D. P.
Tieleman
, and
A. H.
de Vries
,
J. Phys. Chem. B
111
,
7812
(
2007
).
5.
M. F.
Döpke
,
O. A.
Moultos
, and
R.
Hartkamp
,
J. Chem. Phys.
152
,
024501
(
2020
).
6.
A. Z.
Panagiotopoulos
,
Mol. Phys.
61
,
813
(
1987
).
7.
A. Z.
Panagiotopoulos
,
N.
Quirke
,
M.
Stapleton
, and
D. J.
Tildesley
,
Mol. Phys.
63
,
527
(
1988
).
8.
J. L.
Aragones
,
E.
Sanz
, and
C.
Vega
,
J. Chem. Phys.
136
,
244508
(
2012
).
9.
F.
Moučka
,
I.
Nezbeda
, and
W. R.
Smith
,
J. Chem. Phys.
138
,
154102
(
2013
).
10.
I. S.
Joung
and
T. E.
Cheatham
,
J. Phys. Chem. B
113
,
13279
(
2009
).
11.
H.
Jiang
,
I. G.
Economou
, and
A. Z.
Panagiotopoulos
,
Acc. Chem. Res.
50
,
751
(
2017
).
12.
I.
Nezbeda
,
F.
Moučka
, and
W. R.
Smith
,
Mol. Phys.
114
,
1665
(
2016
).
13.
W. R.
Smith
,
I.
Nezbeda
,
J.
Kolafa
, and
F.
Moučka
,
Fluid Phase Equilib.
466
,
19
(
2018
).
14.
C. H.
Pham
,
S. K.
Reddy
,
K.
Chen
,
C.
Knight
, and
F.
Paesani
,
J. Chem. Theory Comput.
13
,
1778
(
2017
).
15.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
5395
(
2013
).
16.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
1599
(
2014
).
17.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
2906
(
2014
).
18.
F. H.
Stillinger
and
A.
Rahman
,
J. Chem. Phys.
60
,
1545
(
1974
).
19.
C.
Vega
and
J. L. F.
Abascal
,
Phys. Chem. Chem. Phys.
13
,
19663
(
2011
).
20.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
21.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
,
J. Chem. Phys.
79
,
926
(
1983
).
22.
I.
Shvab
and
R. J.
Sadus
,
Fluid Phase Equilib.
407
,
7
(
2016
).
23.
P. T.
Kiss
and
A.
Baranyai
,
J. Chem. Phys.
138
,
204507
(
2013
).
24.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
25.
J. R.
Errington
and
A. Z.
Panagiotopoulos
,
J. Phys. Chem. B
102
,
7470
(
1998
).
26.
G.
Lamoureux
,
A. D.
MacKerell
, and
B.
Roux
,
J. Chem. Phys.
119
,
5185
(
2003
).
27.
H. W.
Horn
,
W. C.
Swope
,
J. W.
Pitera
,
J. D.
Madura
,
T. J.
Dick
,
G. L.
Hura
, and
T.
Head-Gordon
,
J. Chem. Phys.
120
,
9665
(
2004
).
28.
R.
Fuentes-Azcatl
and
J.
Alejandre
,
J. Phys. Chem. B
118
,
1263
(
2014
).
29.
R.
Fuentes-Azcatl
,
N.
Mendoza
, and
J.
Alejandre
,
Physica A
420
,
116
(
2015
).
30.
G. A.
Cisneros
,
K. T.
Wikfeldt
,
L.
Ojamäe
,
J.
Lu
,
Y.
Xu
,
H.
Torabifard
,
A. P.
Bartók
,
G.
Csányi
,
V.
Molinero
, and
F.
Paesani
,
Chem. Rev.
116
,
7501
(
2016
).
31.
F.
Moučka
,
I.
Nezbeda
, and
W. R.
Smith
,
J. Chem. Theory Comput.
9
,
5076
(
2013
).
32.
M. B.
Gee
,
N. R.
Cox
,
Y.
Jiao
,
N.
Bentenitis
,
S.
Weerasinghe
, and
P. E.
Smith
,
J. Chem. Theory Comput.
7
,
1369
(
2011
).
33.
I.
Leontyev
and
A.
Stuchebrukhov
,
Phys. Chem. Chem. Phys.
13
,
2613
(
2011
).
34.
Z. R.
Kann
and
J. L.
Skinner
,
J. Chem. Phys.
141
,
104507
(
2014
).
35.
M. P.
Tosi
and
F. G.
Fumi
,
J. Phys. Chem. Solids
25
,
45
(
1964
).
36.
D. E.
Smith
and
L. X.
Dang
,
J. Chem. Phys.
100
,
3757
(
1994
).
37.
G.
Lamoureux
,
E.
Harder
,
I. V.
Vorobyov
,
B.
Roux
, and
A. D.
MacKerell
,
Chem. Phys. Lett.
418
,
245
(
2006
).
38.
I. S.
Joung
and
T. E.
Cheatham
,
J. Phys. Chem. B
112
,
9020
(
2008
).
39.
A. H.
Mao
and
R. V.
Pappu
,
J. Chem. Phys.
137
,
064104
(
2012
).
40.
S.
Reiser
,
S.
Deublein
,
J.
Vrabec
, and
H.
Hasse
,
J. Chem. Phys.
140
,
044504
(
2014
).
41.
P. T.
Kiss
and
A.
Baranyai
,
J. Chem. Phys.
141
,
114501
(
2014
).
42.
R.
Fuentes-Azcatl
and
M. C.
Barbosa
,
J. Phys. Chem. B
120
,
2460
(
2016
).
43.
A. L.
Benavides
,
M. A.
Portillo
,
V. C.
Chamorro
,
J. R.
Espinosa
,
J. L. F.
Abascal
, and
C.
Vega
,
J. Chem. Phys.
147
,
104501
(
2017
).
44.
B.
Widom
,
J. Chem. Phys.
39
,
2808
(
1963
).
45.
P. J.
Lenart
,
A.
Jusufi
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
126
,
044509
(
2007
).
46.
Z.
Abbas
and
E.
Ahlberg
,
J. Solution Chem.
48
,
1222
(
2019
).
47.
E.
Sanz
and
C.
Vega
,
J. Chem. Phys.
126
,
014507
(
2007
).
48.
D.
Frenkel
and
A. J. C.
Ladd
,
J. Chem. Phys.
81
,
3188
(
1984
).
49.
F.
Moučka
,
M.
Lìsal
, and
W. R.
Smith
,
J. Phys. Chem. B
116
,
5468
(
2012
).
50.
F.
Moučka
,
I.
Nezbeda
, and
W. R.
Smith
,
J. Chem. Phys.
139
,
124505
(
2013
).
51.
Z.
Mester
and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
142
,
044507
(
2015
).
52.
H.
Jiang
,
Z.
Mester
,
O. A.
Moultos
,
I. G.
Economou
, and
A. Z.
Panagiotopoulos
,
J. Chem. Theory Comput.
11
,
3802
(
2015
).
53.
Z.
Mester
and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
143
,
044505
(
2015
).
54.
C. W.
Davies
,
J. Chem. Soc.
1938
,
2093
.
55.
M.
Kohns
,
S.
Reiser
,
M.
Horsch
, and
H.
Hasse
,
J. Chem. Phys.
144
,
084112
(
2016
).
56.
W. R.
Smith
,
F.
Moučka
, and
I.
Nezbeda
,
Fluid Phase Equilib.
407
,
76
(
2016
).
57.
M.
Bley
,
M.
Duvail
,
P.
Guilbaud
, and
J.-F.
Dufrêche
,
J. Phys. Chem. B
121
,
9647
(
2017
).
58.
M.
Bley
,
M.
Duvail
,
P.
Guilbaud
, and
J.-F.
Dufrêche
,
J. Phys. Chem. B
122
,
7726
(
2018
).
59.
60.
S. W.
Deleeuw
,
J. W.
Perram
, and
E. R.
Smith
,
Proc. R. Soc. London, Ser. A
373
,
27
(
1980
).
61.
G.
Hummer
,
L. R.
Pratt
, and
A. E.
García
,
J. Phys. Chem.
100
,
1206
(
1996
).
62.
J. P.
Thompson
and
I. C.
Sanchez
,
J. Chem. Phys.
145
,
214103
(
2016
).
63.
J. M.
Young
and
A. Z.
Panagiotopoulos
,
J. Phys. Chem. B
122
,
3330
(
2018
).
64.
M.
Kohns
,
M.
Schappals
,
M.
Horsch
, and
H.
Hasse
,
J. Chem. Eng. Data
61
,
4068
(
2016
).
65.
H.
Jiang
and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
145
,
046101
(
2016
).
66.
F.
Moučka
,
I.
Nezbeda
, and
W. R.
Smith
,
J. Chem. Theory Comput.
11
,
1756
(
2015
).
67.
S.
Mamatkulov
,
M.
Fyta
, and
R. R.
Netz
,
J. Chem. Phys.
138
,
024505
(
2013
).
68.
J. M.
Young
,
C.
Tietz
, and
A. Z.
Panagiotopoulos
,
J. Chem. Eng. Data
65
,
337
(
2020
).
69.
S.
Deublein
,
S.
Reiser
,
J.
Vrabec
, and
H.
Hasse
,
J. Phys. Chem. B
116
,
5448
(
2012
).
70.
T.
Martinek
,
E.
Duboué-Dijon
,
Š.
Timr
,
P. E.
Mason
,
K.
Baxová
,
H. E.
Fischer
,
B.
Schmidt
,
E.
Pluhařová
, and
P.
Jungwirth
,
J. Chem. Phys.
148
,
222813
(
2018
).
71.
L.
Li
,
T.
Totton
, and
D.
Frenkel
,
J. Chem. Phys.
146
,
214110
(
2017
).
72.
L.
Li
,
T.
Totton
, and
D.
Frenkel
,
J. Chem. Phys.
149
,
054102
(
2018
).
73.
J.
Noroozi
and
W. R.
Smith
,
J. Phys. Chem. A
123
,
4074
(
2019
).
74.
J.
Noroozi
and
W. R.
Smith
,
J. Chem. Eng. Data
65
,
1358
(
2020
).
75.
A. L.
Benavides
,
J. L.
Aragones
, and
C.
Vega
,
J. Chem. Phys.
144
,
124504
(
2016
).
76.
J. R.
Espinosa
,
J. M.
Young
,
H.
Jiang
,
D.
Gupta
,
C.
Vega
,
E.
Sanz
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
145
,
154111
(
2016
).
77.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
,
J. Chem. Phys.
144
,
034501
(
2016
).
78.
J.
Kolafa
,
J. Chem. Phys.
145
,
204509
(
2016
).
79.
T.
Yagasaki
,
M.
Matsumoto
, and
H.
Tanaka
,
J. Chem. Theory Comput.
16
,
2460
(
2020
).
80.
I. M.
Zeron
,
J. L. F.
Abascal
, and
C.
Vega
,
J. Chem. Phys.
151
,
134504
(
2019
).
81.
J.
Dočkal
,
M.
Lìsal
, and
F.
Moučka
,
J. Chem. Theory Comput.
16
,
3677
(
2020
).
82.
F.
Moučka
,
J.
Kolafa
,
M.
Lìsal
, and
W. R.
Smith
,
J. Chem. Phys.
148
,
222832
(
2018
).
83.
K. J.
Müller
and
H. G.
Hertz
,
J. Phys. Chem.
100
,
1256
(
1996
).
84.
H. D. B.
Jenkins
and
Y.
Marcus
,
Chem. Rev.
95
,
2695
(
1995
).
85.
S.
Koneshan
,
J. C.
Rasaiah
,
R. M.
Lynden-Bell
, and
S. H.
Lee
,
J. Phys. Chem. B
102
,
4193
(
1998
).
86.
S. H.
Lee
and
J. C.
Rasaiah
,
J. Phys. Chem.
100
,
1420
(
1996
).
87.
A.
Salis
and
B. W.
Ninham
,
Chem. Soc. Rev.
43
,
7358
(
2014
).
88.
G. A.
Orozco
,
O. A.
Moultos
,
H.
Jiang
,
I. G.
Economou
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
141
,
234507
(
2014
).
89.
J. S.
Kim
,
Z.
Wu
,
A. R.
Morrow
,
A.
Yethiraj
, and
A.
Yethiraj
,
J. Phys. Chem. B
116
,
12007
(
2012
).
90.
S.
Yue
and
A. Z.
Panagiotopoulos
,
Mol. Phys.
117
,
3538
(
2019
).
91.
M.
Andreev
,
A.
Chremos
,
J.
de Pablo
, and
J. F.
Douglas
,
J. Phys. Chem. B
121
,
8195
(
2017
).
92.
M.
Andreev
,
J. J.
de Pablo
,
A.
Chremos
, and
J. F.
Douglas
,
J. Phys. Chem. B
122
,
4029
(
2018
).
93.
V.
Agarwal
and
B.
Peters
, “
Solute precipitate nucleation: A review of theory and simulation advances
,” in
Advances in Chemical Physics
, edited by
S. A.
Rice
and
A. R.
Dinner
(
John Wiley & Sons, Inc.
,
2014
), Vol. 155, pp.
97
159
.
94.
D.
Chakraborty
and
G. N.
Patey
,
J. Phys. Chem. Lett.
4
,
573
(
2013
).
95.
G.
Lanaro
and
G. N.
Patey
,
J. Phys. Chem. B
120
,
9076
(
2016
).
96.
G. C.
Sosso
,
J.
Chen
,
S. J.
Cox
,
M.
Fitzner
,
P.
Pedevilla
,
A.
Zen
, and
A.
Michaelides
,
Chem. Rev.
116
,
7078
(
2016
).
97.
G. M.
Torrie
and
J. P.
Valleau
,
J. Comput. Phys.
23
,
187
(
1977
).
98.
A.
Laio
and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
99
,
12562
(
2002
).
99.
C.
Dellago
,
P. G.
Bolhuis
,
F. S.
Csajka
, and
D.
Chandler
,
J. Chem. Phys.
108
,
1964
(
1998
).
100.
R. J.
Allen
,
D.
Frenkel
, and
P. R.
ten Wolde
,
J. Chem. Phys.
124
,
024102
(
2006
).
101.
R. J.
Allen
,
D.
Frenkel
, and
P. R.
ten Wolde
,
J. Chem. Phys.
124
,
194111
(
2006
).
102.
S.
Hussain
and
A.
Haji-Akbari
,
J. Chem. Phys.
152
,
060901
(
2020
).
103.
N. E. R.
Zimmermann
,
B.
Vorselaars
,
D.
Quigley
, and
B.
Peters
,
J. Am. Chem. Soc.
137
,
13352
(
2015
).
104.
N. E. R.
Zimmermann
,
B.
Vorselaars
,
J. R.
Espinosa
,
D.
Quigley
,
W. R.
Smith
,
E.
Sanz
,
C.
Vega
, and
B.
Peters
,
J. Chem. Phys.
148
,
222838
(
2018
).
105.
H.
Jiang
,
A.
Haji-Akbari
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
148
,
044505
(
2018
).
106.
Y.
Gao
,
L. E.
Yu
, and
S. B.
Chen
,
J. Phys. Chem. A
111
,
10660
(
2007
).
107.
H.
Jiang
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
149
,
141102
(
2018
).
108.
D.
Chakraborty
and
G. N.
Patey
,
Chem. Phys. Lett.
587
,
25
(
2013
).
109.
H.
Jiang
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
,
J. Chem. Phys.
150
,
124502
(
2019
).
110.
J.
Anwar
and
P. K.
Boateng
,
J. Am. Chem. Soc.
120
,
9600
(
1998
).
111.
M. N.
Joswiak
,
M. F.
Doherty
, and
B.
Peters
,
Proc. Natl. Acad. Sci. U. S. A.
115
,
656
(
2018
).
112.
M. N.
Joswiak
,
B.
Peters
, and
M. F.
Doherty
,
Cryst. Growth Des.
18
,
6302
(
2018
).
113.
J. R.
Espinosa
,
G. D.
Soria
,
J.
Ramirez
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
,
J. Phys. Chem. Lett.
8
,
4486
(
2017
).
114.
C. J.
Tainter
,
L.
Shi
, and
J. L.
Skinner
,
J. Chem. Theory Comput.
11
,
2268
(
2015
).
115.
H.
Jiang
,
O. A.
Moultos
,
I. G.
Economou
, and
A. Z.
Panagiotopoulos
,
J. Phys. Chem. B
120
,
12358
(
2016
).
116.
R.
Car
and
M.
Parrinello
,
Phys. Rev. Lett.
55
,
2471
(
1985
).
117.
J. Q.
Han
,
L. F.
Zhang
,
R.
Car
, and
E.
Weinan
,
Commun. Comput. Phys.
23
,
629
(
2018
).
118.
J.
Behler
and
M.
Parrinello
,
Phys. Rev. Lett.
98
,
146401
(
2007
).
119.
K. T.
Schütt
,
F.
Arbabzadah
,
S.
Chmiela
,
K. R.
Müller
, and
A.
Tkatchenko
,
Nat. Commun.
8
,
13890
(
2017
).
120.
J.
Sun
,
R. C.
Remsing
,
Y.
Zhang
,
Z.
Sun
,
A.
Ruzsinszky
,
H.
Peng
,
Z.
Yang
,
A.
Paul
,
U.
Waghmare
,
X.
Wu
,
M. L.
Klein
, and
J. P.
Perdew
,
Nat. Chem.
8
,
831
(
2016
).