[1,2]Timothy T. Duignan

1]School of Chemical Engineering, University of Queensland, Brisbane,QLD 4072,
Australia

[2]Queensland Micro- and Nanotechnology Centre, Nathan, Qld 4111, Australia,
Griffith University, Nathan,QLD 4111, Australia

# Accurate, fast and generalisable first principles simulation of aqueous
lithium chloride.

Junji Zhang  Joshua Pagotto  Tim Gould<EMAIL_ADDRESS>[  *

###### Abstract

Electrolyte solutions play a pivotal role throughout chemistry and biology.
For over a century, scientists have therefore sought to accumulate a precise
knowledge and understanding of their thermodynamic, kinetic, and structural
properties. However, the vast majority of electrolyte properties have not _or
cannot_ be determined empirically, despite being key to understanding a huge
range of biological and chemical processes and systems. In this work, we
introduce a long-sought-after solution to this problem and employ it to
develop a detailed understanding of an electrolyte solution: aqueous lithium
chloride. Our solution draws from recent breakthroughs in machine learning,
first principles quantum chemistry and statistical mechanics and lets us
develop and run truly predictive all-atom and coarse-grained simulations using
long-range corrected equivariant neural network potentials (NNP).
Surprisingly, our calculations reveal the formation of Li cation dimers. This
previously unknown species highlights the power of the approach to divulge new
understanding of electrolytes. Key electrolyte properties, including activity
and diffusion coefficients, are determined from first principles and validated
in close agreement with experiment. The training data is a small set (655
frames) of moderate-cost density corrected density functional theory (DC-DFT)
calculations, meaning the approach can be scaled to build a database of
electrolyte solution properties.

###### keywords:

Solutions, liquids, ion-ion interactions, ion diffusivities, thermodynamic
properties, kinetic and structural properties, ab initio molecular dynamics,
equivariant neural network potentials, coarse grained models.

##  1 Introduction

Liquids are the most important phase of matter. Essentially all of biology
occurs in the liquid phase, as does much of chemical engineering and climate
science. They are also notoriously difficult to model or simulate due to their
disordered dynamic nature combined with their sensitivity to complex quantum
mechanical interactions. Recent breakthroughs in machine learning and quantum
chemistry are enabling a paradigm shift towards the prediction of the
properties of liquids starting from first principles.[1, 2, 3, 4]

Electrolyte solutions are a particularly important class of liquids, that are
even more challenging to model due to the inclusion of long range
electrostatic interactions. The prediction of their properties from first
principles is a foundational problem of physical chemistry that has remained
unsolved despite over a century of effort.

Electrolytes are ubiquitous and play an active and central role in a vast
range of important processes and systems. For example, lithium cations are the
primary charge carrier for Li-ion batteries and their chemical equilibria and
diffusivities impact the assembly and performance of these devices.[5, 6]
Remarkably, lithium also exhibits important biochemical effects such as its
efficacy as a treatment for bipolar disorder. However, electrolytes are
important in such a wide range of systems that highlighting individual cases
does not properly convey the scope of applications.

The ability to accurately and efficiently predict the properties of
electrolyte solutions from first principles would be transformative. Firstly,
it would provide direct insight into the molecular scale details that
determine key macroscopic properties. These insights could be used to guide
the design and optimisation of electrolytes for the many important
applications where they play a key role. Secondly, it could be used to
generate large data sets of molecular scale information.[7] Such a database
would be invaluable for training machine learning models that are capable of
generalising to an even wider range of situations,[8] analogous to the way the
protein data bank was leveraged to train Alphafold2.[9]

The foundational model of electrolytes (Debye-Hückel theory) was developed a
century ago this year. It is based on continuum solvent theory and gives a
quantitative description of the thermodynamic and structural properties of
electrolytes but only at millimolar concentrations.[10, 11] At any practically
relevant concentrations, the theory must be corrected using empirically
adjusted parameters. The Pitzer model is the most common example of this
approach and it forms the basis for the calculation of thermodynamic
properties of electrolytes in modern chemical engineering.[12, 13]

This reliance on empirical data to predict the properties of electrolyte
solutions is a critical problem. It means that for any case where good
empirical data does not exist, time consuming experimental characterisation
combined with trial and error is required to understand and optimise
electrolyte solution behaviour. While excellent experimental databases have
been built for some electrolyte solution properties,[13, 14] they are
inherently limited to a relatively small number of cases.[15]

Figure 1: Simulation overview Small, short FPMD simulations are run with CP2K.
Energies/forces are recomputed with the DC-r2SCAN DFA. This data set is then
used to train an equivariant NNP which enables much larger simulations with a
larger cell. Forces and coordinates of the ions alone are output from the NNP-
MD and used to train a coarse grained NNP, which enables much faster MD
simulations and the calculation of key structural, thermodynamic and kinetic
properties.

First Principles Molecular Dynamics (FPMD) based on ab initio or density
functional theory[16, 17, 18, 19] is the only approach that avoids this
dependence on empirical data. However, the extremely large computational
demands of this technique limit its direct applicability to small systems and
a few properties. Additionally, many current density functional approximations
(DFAs) have inaccuracies. [20, 21, 22, 23]

Two exciting new tools enable the transformation of FPMD into a dramatically
more useful tool. The first is density corrected DFT (DC-DFT) which
significantly reduces the errors (delocalisation and self interaction) that
lead standard density functional approximations (DFAs) to inaccurately
describe ions. This method generally uses a Hartree-Fock electron density as
input into the strongly constrained and appropriately normalised (SCAN) DFA.
It has shown particularly promising results recently, including for aqueous
ionic systems. [24, 25, 26, 27, 28, 29]

The second tool is a new generation of Neural Network Potentials (NNP)[30, 31,
32, 33, 34, 35, 36, 37] that incorporates equivariance[38, 39, 40, 41, 42] and
explicit electrostatics.[43, 44, 32, 45, 46] This approach trains a neural
network to reproduce the mapping between coordinates and energies/forces to
massively accelerate FPMD simulation.[30, 31, 32, 33, 34, 35, 36, 37]
Equivaraince means that the rotational symmetries of 3D Euclidean space are
encoded into the neural network. These improvements enable significantly more
accurate, reliable and generalisable NNPs. It has also been demonstrated that
NNPs enable coarse grain molecular dynamics simulation enabling further
acceleration.[47, 48, 49, 50] More generally machine learning potentials are
an incredibly exciting new tool for helping understand many systems, including
electrolytes.[5, 51, 43, 52, 2, 53, 54, 55, 56, 53, 45, 57, 58, 59, 46, 60,
61, 62, 63, 64, 65]

Here, we perform equivariant all atom and coarse grained NNP-MD simulations of
aqueous lithium chloride with explicit long-range electrostatics described by
a continuum solvent model. Density corrected regularised SCAN (DC-r2SCAN) is
used to generate the training data. We observe the remarkable formation of Li
cation dimers and demonstrate excellent agreement with experimental
structural, kinetic, and thermodynamic properties. Figure 1 provides an
overview of the workflow used in this study.

##  2 Results and Discussion

###  2.1 Lithium dimer formation

Figure 2: Lithium dimer formation Lithium oxygen, and lithium lithium RDFs
demonstrate the formation of lithium cation dimers that penetrate into the
first hydration layer and an example snapshot from the MD simulation, showing
the water molecules and nearest chloride ion surrounding the dimer.

The most surprising revelation from the simulations was the formation of
lithium cation dimers, shown in Figure 2, where one lithium penetrates into
the first solvation layer around another lithium ion. The first solvation
layer is defined by the first minima in the Li-O RDF which is 2.8 Å. The
smallest separation of the lithium ions observed in the simulations (2.68 Å)
matches the separation of the neutral covalent dilithium. (2.67 Å) and is much
smaller than the Li-Li distance in LiCl crystal of 3.62 Å. The formation of
this species is particularly surprising given that this ions an archetypal
kosmotrope or water structuring ion. Physically this was believed to
corresponds to the formation of a tightly bonded first solvation layer of
water molecules that was thought to be impenetrable to other ions.[66]

Ironically, the cation dimer appears to be stabilised by lithium’s kosmotropic
behaviour, i.e., the strong structuring of water around the lithium leaves a
region of minimal water density in between the first two hydration layers that
the co-ion can occupy.

This remarkable and counter-intuitive finding may have important implications
for the many biological and chemical systems where lithium plays a critical
role.[6] The transient nature of this pair means that this species could not
feasibly be identified with direct FPMD simulation. CMD simulations show no
indication of it,[67] which is likely attributable to the neglect of charge
transfer and polarisation effects which will significantly mitigate the
electrostatic repulsion. Given the surprising nature of this finding it is
crucial to carefully validate the reliability and robustness of these
simulations as we do in the following sections.

###  2.2 Low training data requirements

The training data for the NNP, trained with NequIP,[39] consists of a
strikingly small training data set of only 655 frames extracted from MD
simulation computed at the DC-r2SCAN level of theory with CP2K.[68, 18, 29]
The 655 frames each contain 4 Li+ and 4 Cl- ions along with 80 water
molecules. The faster r2SCAN level of theory was used to generate the initial
trajectory, which was then re-sampled with DC-r2SCAN.

The total computational cost of generating this dataset was very reasonable,
i.e., on the order of tens of thousands of CPU hours. It is therefore feasbile
to scale this approach to many different electrolyte solutions and build a
large database of properties. By contrast, non-equivariant NNP architectures
generally require much larger training data sets or demand active learning
algorithms to build the data set. Additional details on how the training data
generation are outlined in the Computational Methods section below.

###  2.3 Long, large and accurate simulations

Two NNPs were trained (NNP1 and NNP2) on the DC-r2SCAN data set using
different random seeds for the weight initialisation and for splitting the
data into training and validation sets. NNP-MD was performed with NNP1 on the
original sized cell and 200 frames were extracted from the simulation and
recomputed with DC-r2SCAN level of theory. Figure 3(a) demonstrates that the
NNP is very accurately estimating the forces with an RMSE below 10 meVÅ-1.

Parallel simulations were then run for 2.4 ns each with both NNP1 and NNP2 on
a system six times bigger than the original ( 25.263 Å3 vs. 13.873 Å3)
containing 512 water molecules and 48 ions. Accessing this time and spatial
scale is entirely infeasible with direct FPMD critical, yet is critical for
studying ion-ion correlations which have long range structure that is larger
than the cell size used for FPMD. The NNP-MD generates 200 ps per day for the
larger system on a single V100 GPU.

Figure 3(b),3(c) and 3(d) compares the ion-solvent and ion-ion radial
distribution functions (RDFs) computed with the two separate NNPs showing
excellent agreement and demonstrating convergence and the reproducibility of
the method. The ion-oxygen peak positions of 1.97 and 3.17 Å are in good
agreement with neutron diffraction measurements of 1.96 and 3.15 Å for lithium
and chloride respectively.[69] The first shell coordination number at 4.2 and
6.9 with the NNP-MD are also very consistent with experimental values of 4.1
and 7 from neutron diffraction.[69]

(a) Force correlation plots.

(b) Lithium water RDFs.

(c) Chloride water RDFs.

(d) Ion-ion RDFs.

Figure 3:  Consistency checks of the NNP. (a) Comparison of the forces
computed at the DC-DFT level of theory with the predictions of the NNP for 200
frames extracted from a NNP-MD simulation. (b),(c) and(d) compare the RDFs
predictions of two NNPs trained with different initial seeds demonstrating
consistency. Good agreement with neutron diffraction measurements of the peak
in the ion oxygen RDF is also observed.

###  2.4 Generalisability to lower concentrations

It is important to validate the NNPs description of the water-water
interactions. The structural properties of water, such as the oxygen-oxygen
RDF, computed at 2.5 M are difficult to verify experimentally. However, to
test the NNPs description of the pure water interactions we can run a pure
water simulation. Remarkably Figure 4(a) shows essentially perfect agreement
with experimental X-ray diffraction measurements. The small difference in the
first peak height is attributable to the neglect of quantum nuclear effects,
which is verified by comparison with MB-Pol, a state of the art water model.

(a) Oxygen-oxygen RDF.

(b) Ion-ion RDFS

(c) Ion-ion RDFS 2.5 M

(d) Ion-ion RDFS 1.3 M

Figure 4: (a) Comparison of the predicted O-O RDF with a pure water NNP-MD
simualtion with experimental X-Ray Diffraction measurements[70] and the MB-pol
water model.[71] (b) Comparison of ion-ion RDFs at the concentration of 2.5 M
with a lower concentration (1.3 M) which the model was not trained on. (c,d)
Comparison of the all atom and coarse grained ion-ion RDFs predicted at 2.5 M
(c) and 1.3 M (d) showing good agreement.

The fact that without any explicit effort to fit to or reproduce the
properties of pure water we arrive at a model that gives such accurate
structural predictions is particularly promising and demonstrates the
generalisability of this approach.

To further test the generalisability of the method simulations were performed
at 1.3 M using the same NNPs. Similarly stable simulations were observed with
no physical inaccuracies. Most promisingly an increasing Debye-Hückel
screening length is visible in these simulations as the concentration
decreases. The fact that NNP-MD can reproduce this effect despite using
training data from a single concentration is promising.

###  2.5 Coarse Graining

The all atom NNP-MD is much faster than FPMD and capable of simulating
experimentally relevant timescales but the computational cost is still non-
trivial and more expensive than most CMD approaches. To remedy this we build a
coarse grained model of the electrolyte solutions, more specifically we
integrate out the solvent degrees of freedom resulting in a continuum or
implicit solvent model. To do this we use the NNP to learn the potential of
mean force, which is a free energy surface. The problem of learning a free
energy surface is very similar to learning the full all-atom potential energy
surface and so also benefits from equivariance as has recently been
demonstrated.[48] As before we compute the long range electrostatic
interactions separately using the known analytical expression.

In practice to do this we extract the coordinates and forces for the ions
alone for 24,000 frames extracted from a 2.4 ns NNP-MD all atom simulation. A
larger data set is required to sufficiently converge the averaging over the
solvent degrees of freedom. It would not be feasible to generate such a large
training data set with FPMD directly but is straightforward with NNP-MD. The
coarse grain NNP requires many fewer weights and trains very quickly in
comparision to the all atom NNP due to the much simpler energy surface.

The coarse grained MD can reproduce the RDFs from the all atom MD accurately
as shown in Figure 4(c). It is also reasonably generalisable as shown in
Figure 4(d) where it reproduces the RDF at 1.3 M including the decrease in
long range screening. This is impressive considering that we only trained on
data from a single concentration.

The coarse grained NNP-MD is orders of magnitude faster than the all atom NNP
and requires trivial computational resources to fully converge, i.e, tens of
CPU hours. This hierarchical layering of NNPs where a coarse grained NNP is
trained on a lower level data from all atom NNP demonstrates a promising
general approach to the long standing challenge of connecting scales in
molecular simulation. [72, 73]

##  3 Empirical validation

While there is no way to directly validate the ion-ion RDFs, as they cannot be
accurately determined experimentally, their reliability can be indirectly
validated by comparison with thermodynamic properties. Specifically, the
activity coefficients as a function of concentration have been shown to be
highly sensitive to the strength of the ion-ion interactions. The accurate
prediction of these quantities gave the original validation of Debye-Hückel
theory. The ability to determine these properties accurately is also of
immense practical importance for chemical engineering.

Here, we use Kirkwoof-Buff theory to compute the derivative of the activities
at 2.5 and 1.3 M, with good experimental agreement as shown in Figure 5(a) for
both the all atom MD and coarse grained MD RDFs. The coarse grained
predictions have negligible error bars due to their better convergence. The
all atom MD is used for the ion-water RDFs, which converge much more quickly
than the ion-ion RDFs and make a very small contribution to the total.

Diffusivites of Li+ and Cl- ions as well as water molecules are also in good
agreement with experiment as shown in Figure 5(b). This is particularly
impressive as kinetic properties depend on accurately assessing barrier
heights in the potential energy surface. These will be poorly represented in
the training data set as it is extracted from equilibrium molecular dynamics.

A direct relative of the activity coefficients, the osmotic coefficients, can
also be approximately computed using a modification of Debye-Hückel’s theory.
The key requirement is the infinite dilution potential of mean force (PMF),
also known as the pair potential. This is a fundamental quantity that can be
inserted into the Poisson-Boltzmann equation and solved to compute
electrostatic potential and RDFs at any concentration with trivial
computational cost. A virial expression can then be used to convert these to
osmotic coefficients.

Computing the infinite dilution PMF is normally an extremely challenging task
requiring enhanced sampling techniques, such as umbrella sampling with very
low concentration simulations. Here we can trivially extract this quantity
from the coarse grained NNP by simply computing the interaction energy between
the three ion pairs in isolation. These infinite dilution PMFs are shown in
Figure 5(c) with the long range Coulomb term removed. A correction term is
included to account for the variation in the dielectric constant as a function
of concentration between the training concentration of 2.5 M and infinite
dilution.

(a) Activity coefficient derivatives.

(b) Diffusivities

(c) Short range infinite dilution PMF.

(d) Osmotic coefficients.

Figure 5: (a) Comparison of activity coefficient derivatives predicted with
Kirkwood-Buff theory using both all atom and coarse grained MD at two
concentrations. (b) Comparison of the computed diffusivities with experimental
values[74] for water and both ions at two concentrations and infinite dilution
also showing good experimental agreement. (c) The infinite dilution potential
of mean force between the two ions computed with the coarse grained NNP. (d)
Comparison of osmotic coefficients computed with the short ranged infinite
dilution PMF and the modified Poisson-Boltzmann equation.

Figure 5(d) shows the theoretically predicted osmotic coefficients
demonstrating good agreement with experiment. These can be directly converted
to activity coefficients using the Gibbs-Duhem relation. The experimental
values are converted from Lewis–Randall to Macmillan–Mayer scheme as
required.[75, 76]

##  4 Conclusion and future work.

In summary, we have demonstrated the accurate prediction of the structural,
kinetic and thermodynamic properties of aqueous LiCl electrolyte solution. Our
approach combines the exciting new tools of equivariant NNPs and DC-DFT with
the century old mathematical theory of electrolytes. Despite its remarkable
predicitive ability, it has very low training data requirements and requires
relatively few computational resources.

The prediction of these properties, particularly activity coefficients, has
been a key goal of physical chemistry for over a century. This is firstly
because they play a critical role in huge range of important applications and
secondly because the reproduction of these properties is an important
validation of the reliability and accuracy of these methods, which is a
necessary prerequisite to the accurate simulation of many more complex
systems.

The next step is to begin scaling this method to build a database of
properties of electrolytes across a much wider range of conditions and
compositions than exists currently. A key focus should be on important
electrolytes that are particularly difficult to characterise experimentally
such as pure lithium bicarbonate which immediately speciates into a mixture of
carbonates in solution, meaning even its most basic properties have never been
measured. More advanced techniques for generating the data sets such as meta-
dynamics[77] and active learning[61, 65] will likely be required for these
more complex systems.

Besides being of direct practical use, this database can also provide the
training data for a new generation of AI models analogous to Alphafold2’s use
of the PDB.[9] This should enable the rapid prediction of the properties of
electrolyte solutions, without requiring direct simulation. Such a database
could also be used to train generative models for designing new
electrolytes.[78]

##  5 Computational Methods

###  5.1 FPMD (CP2K)

We used Born-Oppenheimer ab initio molecular dynamics simulations within a
constant volume NVT (300 K) ensemble with periodic boundary conditions. The
CP2K simulation suite (http:www.cp2k.org) containing the _QuickStep_ module
for the DFT calculations [68, 18] was used with a 0.5 fs time step. We used a
double $\zeta$ basis set that has been optimized for the condensed phase[79]
in conjunction with GTH pseudopotentials [80] optimised for SCAN[81, 82] and a
1200 Ry cutoff.[83, 84] A Nosé-Hoover thermostat was attached to every degree
of freedom to ensure equilibration. [85]

An $\approx$ 10 ps simulation was run consisiting of 4 lithium ions, 4
chloride ions and 80 water molecules in a cell of fixed dimensions of 13.873
Å3, corresponding to an electrolyte concentration of 2.5 M. The cell size was
adjusted to match the experimental density at this concentration. The initial
simulation used the r2SCAN DFA.[24]

458 frames were extracted from this initial simulation in addition to 197
frames from a 1 ns NNP-MD simulation trained on it. (see below for details)
The forces were then recomputed on these 655 frames using the density
corrected r2SCAN DFA. [27, 26] This method has recently been implemented in
CP2K[29]. The auxiliary density matrix method (ADMM) was employed to improve
the scaling of the four-center two-electron integrals.[86] The Schwarz
integral screening threshold was set to 10-5 units. A contracted auxiliary
basis set (cFIT3) was used to construct the auxiliary density matrix.

###  5.2 NNP fitting (NequIP)

The training data set was generated as follows: 1) 458 (70%) frames were
sampled directly from a 10 ps FPMD run using r2SCAN; 2) a trial NNP was fit to
the r2SCAN forces, and used to generate 1 ns of NNP-MD data with the FPMD box
size; 3) the remaining 197 (30%) frames were then sampled from the 1 ns NNP-MD
run, to improve the diversity of training data for the more sophisticated NNP
resulting in total of 655 frames.

The trial r2SCAN was fit to 985 frames. The 655 frames of DC-r2SCAN
calculations from both the FPMD and the NNP-MD were used to train two new
neural networks, referred to as NNP1 and NNP2 above.

As our NNP only has access to local information, (short range, $<5$ Å)[44] the
long range electrostatic ion-ion interactions were removed from the forces and
energies prior to training. These were computed using a dielectrically
screened Coulomb interaction.[87] They were then added back in during all NNP-
MD simulations.

The same hyperparameters were used for all of the all atom NNPs. More
specifically 100:1 weighting on forces vs energies was used in the default
loss function.[39] We decrease the initial learning rate of 0.01 by a decay
factor of 0.5 whenever the validation RMSE in the forces has not seen an
improvement for five epochs. Training was stopped when the learning rate is
smaller than 10-5. The model with the lowest validation error was used for
simulations. A radial cutoff distance of 5 Å was used. Three layers of
interaction blocks were used with the maximum $l$ set to 2 each with 16
features. Only even parity was used. Invariant neurons for the radial netowrk
was set to 32. All the other parameters were set to the defaults. A 80:20
training validation split was used throughout.

For all the NNPs the Coulomb interactions screened by the dielectric constant
of water (78.3) between all the ions were subtracted before training. These
were calculated with LAMMPS by placing appropriately screened charges on the
ions to reproduce dielectric screening of 78.3 and with the particle-particle
particle-mesh method.[88]

To train the NNP for the coarse grain MD we found that reducing the number of
parameters was important to provide stability. Two layers of interaction
blocks were used with the maximum $l$ set to one each with 8 features. Only
even parity was used. Invariant neurons for the radial network was set to 8.
All the other parameters were set to the defaults. The radial cut off was also
extended to 10 Å to provide longer range interactions.

Figure. 6(a) and 6(b) show the learning curves. RMSE on the validation set for
the first all atom NNP (NNP1) was 8.9 meV/Å for the forces and 0.11 meV for
the energies. For NNP2 it was 9.7 meV/Å for the forces and 0.12 meV for the
energies. The RMSE on the ions is much higher with the coarse grained model as
expected due to the neglect of the solvent. They were 299 meV/Å for the forces
and 38 meV for the energies. Remarkably the error on initialisation was 302
meV/Å and 91 meV respectively meaning that training of the NNP only removed 3
meV/Å in error and yet this is enough to reliably reproduce the ion-ion RDFs.

(a) All atom NNP

(b) Coarse grained NNP

Figure 6: Learning curves

###  5.3 NNP MD (NequIP/LAMMPS)

The NequIP plugin for LAMMPS[89] was used to perform several NVT simulations
at 300 K for 2.4 ns. A Nosé-Hoover thermostat was attached to every degree of
freedom to ensure equilibration [85]. The long-range Coulomb interactions were
added to the simulation using LAMMPS hybrid overlay method. No initial data
was discarded as the initial frame was taken from the end of the AIMD
simulation.

We use constant volume simulations with the experimental density to estimate
the cell size to avoid relying on experimental information. The 2.5 M
simulations in the larger cell size had dimensions of 25.263 Å3 and contained
48 ions and 512 water molecules. The 1.3 M simulations had a 25.053 Å3 cell
size with 24 ions and 512 water molecules. Three 2.4 ns simulations were run
at 2.5 M and four at 1.2 M to ensure convergence. Most simulations were very
stable with no noticeable non-physical behaviour except for one simulation at
2.5 M where water molecule deprotonated with the hydrogen forming a tight bond
with lithium. This trajectory was not used.

VMD[90] was used to create the RDFs and images in Fig. 1.

The infinite dilution PMFs were computed in LAMMPS by simply computing the
total energy of the system of two ions in a 25.263 Å3 as a function of
distance. At short distances, where the NNP has no data as the ions do not
approach closely, the NNP can oscillate randomly. At these points an
increasing extrapolation was used to ensure that the infinite dilution PMF
didn’t go negative again.

###  5.4 Modified Poisson-Boltzmann Equation (MPBE)

The dielectric Coulomb potential was removed analytically to determine the
short range contribution to the infinite diultion PMF and a correction was
included in the short range potential, i.e., between 2 and 10 Å to account for
the change in dielectric constant between 2.5 M and infinite dilution. The
correction is given by:

$\frac{qe}{4\pi\epsilon_{\text{w}}r}-\frac{qe}{4\pi\epsilon_{\text{s}}r}$ (1)

where the dielectric constant of bulk water was given by:
$\epsilon_{\text{w}}=78.3$ and the dielectric constant of the solution at 2.5
M is given by $\epsilon_{\text{s}}=\frac{78.3}{1+\beta\rho}$. This is a
standard expression used to compute the concentration dependence of the
dielectric constant.[91] We use the value of $\beta=0.19$.

The infinite dilution PMF was then used in the MBPE to compute the
electrostatic potential and RDFs at 17 concentrations between 0.01 and 4 M.
The MPBE is given by:

$-\epsilon_{r}\epsilon_{o}\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\phi(r)}{dr}\right)=\sum_{i}q_{i}\rho_{i}(r)$
(2)

Where the density is given by:

$\rho_{i}(r)=\rho_{i}\exp{\left[-\beta\left(q_{i}\phi(r)+W_{\text{SR}}(r)\right)\right]}$
(3)

$\frac{\rho_{i}(r)}{\rho_{i}}$ corresponds to the RDF. The boundary conditions
used to solve the MPBE are the electric field of a point charge at 2 Å and the
electric potential of 0 at large separations (60 Å). These settings ensure
reliable numerical solutions are found using the shooting method implemented
in Mathematica.[92] Smoothing was applied to the PMFs to improve the
convergence of the numerical differential equation solution.

Note that there are two ways to solve the MPBE one with the cation as the
central ion and one with the anion. The choice does not materially effect the
predicted osmotic coefficients though.

###  5.5 Osmotic coefficients calculation

The MPBE can determine RDFs down to infinite dilution with minimal
computational demands. $W_{\text{SR}}(r)$ can then be input into the virial
expression for the osmotic coefficients to estimate them at many
concentrations.[93, 67, 94]

$\phi(\rho)=1-\frac{\pi}{3}\rho\sum_{i,j}\int_{0}^{\infty}g_{ij}(r)\frac{dW_{ij}(r)}{dr}r^{3}dr$
(4)

where $W_{ij}(r)$ is the ion-ion infinite dilution PMF, i.e.,
$W_{\text{SR}}(r)$ plus the Coulomb interaction, i.e., the full infinite
dilution PMF. This gives reasonable agreement with osmotic coefficients as
shown in the main text.

Experimental densities[95, 96, 97] were used to convert theoretical osmotic
coefficients from the McMillan-Mayer to the Lewis-Randal level to enable the
correct thermodynamic comparison.[75, 76] Experimental values are obtained
from the Pitzer equations and converted to molarity again using experimental
densities.[98] A fixed dielectric constant was used in the MPBE as this
corresponds best with underlying McMillan-Mayer level the theory is based on,
i.e., fixed solvent chemical potential.

The Gibbs-Duhem relationship between osmotic and activity coefficients could
be used to convert osmotic coefficients into activities if desired.

###  5.6 Kirkwood-Buff theory calculations.

Kirkwood-Buff theory [99] was used to compute the derivatives of the
activities from the RDFs with the following expressions:

$\frac{1}{1+\rho\left(G_{\text{cc}}-G_{\text{co}}\right)}$ (5)

where $\rho$ is the ion density. $G_{\text{cc}}$ is given for a monovalent
ions by:

$G_{\text{cc}}=\frac{1}{4}\left(G_{++}+G_{--}+2G_{+-}\right)$ (6)

and $G_{\text{co}}$ is given by:

$G_{\text{co}}=\frac{1}{2}\left((G_{+\text{O}}+G_{-\text{O}})\right)$ (7)

where the $G$ refers to the Kirkwood-Buff integrals:

$G_{ij}=\int_{0}^{\infty}\left(g_{ij}(r)-1\right)r^{2}\,dr$ (8)

The integrals were cutoff at half the box size and the RDFs were normalised by
their average value around 1 Å of the cutoff to ensure they went to 1. For the
coarse grained MD results the ion-solvent integrals were computed using the
coarse-grained MD.

###  5.7 Diffusion coefficients calculation

Diffusion coefficients were computed from the mean squared displacements (MSD)
of the water molecule and lithium and chloride ions in our NNP MD
trajectories. This conversion was carried out using the diffusion coefficient-
MSD relationship described below:

$D=\frac{\text{MSD}}{6t}$ (9)

The results were finally adjusted by finite size corrections [100]. Here, we
have used the experimental value (0.888 mPas) for the viscosity of pure water
when determining the finite size correction. Experimental values were obtained
from Ref. 74.

##  6 Acknowledgements

We acknowledge Stefan Vuckovic, Debra Bernhardt, Christopher Mundy, Gregory
Schenter, Simon Batzner, Alby Musaelian, Sophie Baker, Alister Page, Juerg
Hutter and Andrey Bliznyuk for helpful discussions and computational
assistance. TTD was supported by an Australian Research Council (ARC)
Discovery Project (DP200102573) and DECRA Fellowship (DE200100794). TG was
supported by an Australian Research Council (ARC) Discovery Project
(DP200100033) and Future Fellowship (FT210100663). This research was
undertaken with the assistance of resources and services from the National
Computational Infrastructure (NCI), which is supported by the Australian
Government and with the assistance of resources from QCIF.
(http://www.qcif.edu.au)

## References

  *   * [1] Deringer, V. L. _et al._ Origins of structural and electronic transitions in disordered silicon.  _Nat_ 589, 59–64 (2021).
  * [2] Galib, M. & Limmer, D. T.  Reactive uptake of N2O5 by atmospheric aerosol is dominated by interfacial processes.  _Science (New York, N.Y.)_ 371, 921–925 (2021).
  * [3] Kapil, V. _et al._ The first-principles phase diagram of monolayer nanoconfined water.  _Nat_ 609, 512–516 (2022).
  * [4] Bore, S. L. & Paesani, F.  Realistic phase diagram of water from “first principles” data-driven quantum simulations.  _Nat. Commun._ 14, 3349 (2023).
  * [5] Hellström, M. & Behler, J.  Structure of aqueous NaOH solutions: Insights from neural-network-based molecular dynamics simulations.  _Phys. Chem. Chem. Phys._ 19, 82–96 (2017).
  * [6] Chen, X. & Zhang, Q.  Atomic Insights into the Fundamental Interactions in Lithium Battery Electrolytes.  _Accounts Chem. Res._ 53, 1992–2002 (2020).
  * [7] Gregory, K. P., Elliott, G. R., Wanless, E. J., Webber, G. B. & Page, A. J.  A quantum chemical molecular dynamics repository of solvated ions.  _Sci. Data_ 9, 430 (2022).
  * [8] Winter, B., Winter, C., Esper, T., Schilling, J. & Bardow, A.  SPT-NRTL: A physics-guided machine learning model to predict thermodynamically consistent activity coefficients.  _Fluid Phase Equilibr._ 568, 113731 (2023).
  * [9] Jumper, J. _et al._ Highly accurate protein structure prediction with AlphaFold.  _Nat_ 596, 583–589 (2021).
  * [10] Debye, P. & Hückel, E.  Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und Verwandte Erscheinungen [The Theory of Electrolytes. I. Lowering of Freezing Point and Related Phenomena].  _Phys. Z._ 24, 185–206 (1923).
  * [11] Robinson, R. A. & Stokes, R. H.  _Electrolyte Solutions_ (Butterworth & Co., Devon, 1959).
  * [12] Pitzer, K. S.  Thermodynamics of electrolytes. I. Theoretical basis and general equations.  _J. Phys. Chem._ 77, 268–277 (1973).
  * [13] May, P. M. & Rowland, D.  Thermodynamic Modeling of Aqueous Electrolyte Systems: Current Status.  _J. Chem. Eng. Data_ 62, 2481–2495 (2017).
  * [14] Rowland, D. & May, P. M.  Progress in Aqueous Solution Modelling: Better Data and Better Interfaces.  _J. Solution Chem._ 48, 1066–1078 (2019).
  * [15] Vaque Aura, S. _et al._ Data Analysis for Electrolyte Systems: A Method Illustrated on Alkali Halides in Water.  _J. Chem. Eng. Data_ 66, 2976–2990 (2021).
  * [16] Marx, D. & Hutter, J.  _Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods_ (Cambridge University Press, Cambridge, 2009).
  * [17] Baer, M. D. & Mundy, C. J.  Toward an Understanding of the Specific Ion Effect Using Density Functional Theory.  _J. Phys. Chem. Lett._ 2, 1088–1093 (2011).
  * [18] Kühne, T. D. _et al._ CP2K: An Electronic Structure and Molecular Dynamics Software Package I. Quickstep: Efficient and Accurate Electronic Structure Calculations.  _J. Chem. Phys._ 152, 1–51 (2020).
  * [19] Duignan, T. T., Kathmann, S. M., Schenter, G. K. & Mundy, C. J.  Toward a First-Principles Framework for Predicting Collective Properties of Electrolytes.  _Accounts Chem. Res._ 54, 2833–2843 (2021).
  * [20] Riera, M., Mardirossian, N., Bajaj, P., Götz, A. W. & Paesani, F.  Toward chemical accuracy in the description of ion–water interactions through many-body representations. Alkali-water dimer potential energy surfaces.  _J. Chem. Phys._ 147, 161715 (2017).
  * [21] Paesani, F., Bajaj, P. & Riera, M.  Chemical accuracy in modeling halide ion hydration from many-body representations.  _Advances in Physics: X_ 4, 1631212 (2019).
  * [22] Wagle, K. _et al._ Self-Interaction Correction in Water – Ion Clusters.  _J. Chem. Phys._ 154, 094302 (2021).
  * [23] Duignan, T. T., Schenter, G., Mundy, C. J. & Zhao, X. S.  A method for accurately predicting solvation structure.  _J. Chem. Theory Comput._ 16, 5401–5409 (2020).
  * [24] Furness, J. W., Kaplan, A. D., Ning, J., Perdew, J. P. & Sun, J.  Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation.  _J. Phys. Chem. Lett._ 11, 8208–8215 (2020).
  * [25] Dasgupta, S., Lambros, E., Perdew, J. P. & Paesani, F.  Elevating Density Functional Theory to Chemical Accuracy for Water Simulations through a Density-Corrected Many-Body Formalism.  _Nature comm_ 12, 6359 (2021).
  * [26] Palos, E., Caruso, A. & Paesani, F.  Consistent Density Functional Theory-Based Description of Ion Hydration Through Density-Corrected Many-Body Representations.  _ChemRXiv_ 10.26434/chemrxiv–2023–vp6ns–v2 (2023).
  * [27] Sim, E., Song, S., Vuckovic, S. & Burke, K.  Improving Results by Improving Densities: Density-Corrected Density Functional Theory.  _J. Am. Chem. Soc._ 144, 6625–6639 (2022).
  * [28] Song, S. _et al._ Extending density functional theory with near chemical accuracy beyond pure water.  _Nature Communications_ 14, 799 (2023).
  * [29] Belleflamme, F. & Hutter, J.  Radicals in aqueous solution: Assessment of density-corrected SCAN functional.  _Phys. Chem. Chem. Phys._ 10.1039.D3CP02517A (2023).
  * [30] Mater, A. C. & Coote, M. L.  Deep Learning in Chemistry.  _J. Chem. Inf. Model._ 59, 2545–2559 (2019).
  * [31] Noé, F., Tkatchenko, A., Müller, K.-R. & Clementi, C.  Machine Learning for Molecular Simulation.  _Annu. Rev. Phys. Chem._ 71, 361–390 (2020).
  * [32] Behler, J.  Four Generations of High-Dimensional Neural Network Potentials.  _Chem. Rev._ 121, 10037–10072 (2021).
  * [33] White, A. D.  Deep Learning for Molecules and Materials.  _Living Journal of Computational Molecular Science_ 3, 1499 (2022).
  * [34] Kocer, E., Ko, T. W. & Behler, J.  Neural Network Potentials: A Concise Overview of Methods.  _Annu. Rev. Phys. Chem._ 73, 163–186 (2022).
  * [35] Wen, T., Zhang, L., Wang, H., E, W. & Srolovitz, D. J.  Deep potentials for materials science.  _Materials Futures_ 1, 022601 (2022).
  * [36] Kolluru, A. _et al._ Open Challenges in Developing Generalizable Large Scale Machine Learning Models for Catalyst Discovery.  _arXiv_ arXiv:2206.02005 (2022).
  * [37] Tokita, A. M. & Behler, J.  Tutorial: How to Train a Neural Network Potential.  _arXiv_ arXiv:2308.08859 (2023).
  * [38] Satorras, V. G., Hoogeboom, E. & Welling, M.  E(n) Equivariant Graph Neural Networks.  _arXiv_ arXiv:2102.09844 (2021).
  * [39] Batzner, S. _et al._ E(3)-Equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials.  _Nat. Commun._ 13, 2453 (2022).
  * [40] Batatia, I. _et al._ The Design Space of E(3)-Equivariant Atom-Centered Interatomic Potentials.  _arXiv_ arXiv:2205.06643 (2022).
  * [41] Tholke, P. & Fabritiis, G. D.  TorchMD-NET: Equivariant Transformers for Neural Network based Molecular Potentials.  _arXiv_ arXiv:2202.02541 (2022).
  * [42] Liao, Y.-l. & Smidt, T.  Equiformer: Equivariant Graph Attention Transformer for 3D Atomistic Graphs.  _ArXiV_ arXiv:2206.11990 (2022).
  * [43] Yao, K., Herr, J. E., Toth, D. W., Mckintyre, R. & Parkhill, J.  The TensorMol-0.1 model chemistry: A neural network augmented with long-range physics.  _Chem. Sci._ 9, 2261–2269 (2018).
  * [44] Yue, S. _et al._ When do short-range atomistic machine-learning models fall short?  _J. Chem. Phys._ 154, 034111 (2021).
  * [45] Gao, A. & Remsing, R. C.  Self-Consistent Determination of Long-Range Electrostatics in Neural Network Potentials.  _Nat. Commun._ 13, 1572 (2021).
  * [46] Jacobson, L. D. _et al._ Transferable Neural Network Potential Energy Surfaces for Closed-Shell Organic Molecules: Extension to Ions.  _J. Chem. Theory Comput._ 18, 2354–2366 (2022).
  * [47] Husic, B. E. _et al._ Coarse graining molecular dynamics with graph neural networks.  _J. Chem. Phys._ 153, 194101 (2020).
  * [48] Loose, T. D., Sahrmann, P. G., Qu, T. S. & Voth, G. A.  Coarse-Graining with Equivariant Neural Networks: A Path Towards Accurate and Data-Efficient Models.  _arXiv_ arxiv:2309.09290 (2023).
  * [49] Majewski, M. _et al._ Machine learning coarse-grained potentials of protein thermodynamics.  _Nat. Commun._ 14, 5739 (2023).
  * [50] Wilson, M. O. & Huang, D. M.  Anisotropic molecular coarse-graining by force and torque matching with neural networks.  _J. Chem. Phys._ 159, 024110 (2023).
  * [51] Smith, J. S., Isayev, O. & Roitberg, A. E.  ANI-1: An extensible neural network potential with DFT accuracy at force field computational cost.  _Chem. Sci._ 8, 3192–3203 (2017).
  * [52] Onat, B., Cubuk, E. D., Malone, B. D. & Kaxiras, E.  Implanted neural network potentials: Application to Li-Si alloys.  _Phys. Rev. B_ 97, 1–9 (2018).
  * [53] Zhang, L., Wang, H., Car, R. & Weinan, E.  Phase Diagram of a Deep Potential Water Model.  _Phys. Rev. Lett._ 126, 236001 (2021).
  * [54] Shi, Y., Doyle, C. C. & Beck, T. L.  Condensed Phase Water Molecular Multipole Moments from Deep Neural Network Models Trained on Ab Initio Simulation Data.  _J. Phys. Chem. Lett._ 12, 10310–10317 (2021).
  * [55] Guo, Z. _et al._ Extending the limit of molecular dynamics with ab initio accuracy to 10 billion atoms.  _arXiv_ arXiv:2201.01446 (2022).
  * [56] Zhang, C., Yue, S., Panagiotopoulos, A. Z., Klein, M. L. & Wu, X.  Dissolving salt is not equivalent to applying a pressure on water.  _Nat. Commun._ 13, 822 (2022).
  * [57] Zhang, J., Pagotto, J. & Duignan, T. T.  Towards predictive design of electrolyte solutions by accelerating _ab initio_ simulation with neural networks.  _J. Mater. Chem.A_ 10, 19560–19571 (2022).
  * [58] Malosso, C., Zhang, L., Car, R., Baroni, S. & Tisi, D.  Viscosity in water from first-principles and deep-neural-network simulations.  _npj Comput. Mater._ 8, 139 (2022).
  * [59] Shi, Y., Lam, S. T. & Beck, T. L.  Quasichemical theory and quantum simulation based deep-learning potentials for modeling molten salt thermodynamics.  _Chem. Sci._ 13, 8265–8273 (2022).
  * [60] Dajnowicz, S. _et al._ High-Dimensional Neural Network Potential for Liquid Electrolyte Simulations.  _The Journal of Physical Chemistry B_ 126, 6271–6280 (2022).
  * [61] Wang, F. & Cheng, J.  Automated Workflow for Computation of Redox Potentials, Acidity Constants and Solvation Free Energies Accelerated by Machine Learning.  _J. Chem. Phys._ 157, 024103 (2022).
  * [62] Zhang, C., Yue, S., Panagiotopoulos, A. Z., Klein, M. L. & Wu, X.  Why Dissolving Salt in Water Decreases Its Dielectric Permittivity.  _Phys. Rev. Lett._ 131, 076801 (2023).
  * [63] Avula, N. V. S., Klein, M. L. & Balasubramanian, S.  Understanding the Anomalous Diffusion of Water in Aqueous Electrolytes Using Machine Learned Potentials.  _arXiv_ arXiv:2307.15576 (2023).
  * [64] Jacobson, L., Stevenson, J., Ramezanghorbani, F., Dajnowicz, S. & Leswing, K.  Leveraging Multitask Learning to Improve the Transferability of Machine Learned Force Fields.  _ChemRXiv_ 10.26434/chemrxiv–2023–8n737 (2023).
  * [65] Guo, J. _et al._ AL4GAP: Active Learning Workflow for generating DFT-SCAN Accurate Machine-Learning Potentials for Combinatorial Molten Salt Mixtures.  _ChemRXiv_ 10.26434/chemrxiv–2023–wzv3q (2023).
  * [66] Collins, K. D.  Charge Density-Dependent Strength of Hydration and Biological Structure.  _Biophys. J._ 72, 65–76 (1997).
  * [67] Kalcher, I. & Dzubiella, J.  Structure-Thermodynamics Relation of Electrolyte Solutions.  _J. Chem. Phys._ 130, 134507 (2009).
  * [68] VandeVondele, J. _et al._ Quickstep: Fast and accurate density functional calculations using a mixed gaussian and plane waves approach.  _Comput. Phys. Commun._ 167, 103–128 (2005).
  * [69] Marcus, Y.  Effect of ions on the structure of water: Structure making and breaking.  _Chem. Rev._ 109, 1346–1370 (2009).
  * [70] Skinner, L. B. _et al._ Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range.  _J. Chem. Phys._ 138, 074506 (2013).
  * [71] Medders, G. R., Babin, V. & Paesani, F.  Development of a “first-principles” water potential with flexible monomers. III. Liquid phase properties.  _J. Chem. Theory Comput._ 10, 2906–2910 (2014).
  * [72] Wang, J. _et al._ Machine Learning of Coarse-Grained Molecular Dynamics Force Fields.  _ACS Cent. Sci._ 5, 755–767 (2019).
  * [73] E, W., Lei, H., Xie, P. & Zhang, L.  Machine learning-assisted multi-scale modeling.  _J. Math. Phys._ 64, 071101 (2023).
  * [74] Tanaka, K. & Nomura, M.  Measurements of tracer diffusion coefficients of lithium ions, chloride ions and water in aqueous lithium chloride solutions.  _Journal of the Chemical Society, Faraday Transactions 1: Physical Chemistry in Condensed Phases_ 83, 1779 (1987).
  * [75] Friedman, H. L.  Lewis-Randall to McMillan-Mayer Conversion for the Thermodynamic Excess Functions of Solutions. Part II. Excess Energy and Volume.  _J. Solution Chem._ 1, 413–417 (1972).
  * [76] Simonin, J.-P.  Study of Experimental-to-McMiIlan-Mayer Conversion of Thermodynamic Excess Functions.  _Journal of the Chemical Society, Faraday Transactions_ 92, 3519–3523 (1996).
  * [77] Herr, J. E., Yao, K., McIntyre, R., Toth, D. W. & Parkhill, J.  Metadynamics for training neural network model chemistries: A competitive assessment.  _J. Chem. Phys._ 148, 241710 (2018).
  * [78] Kaufman, B. _et al._ COATI: Multi-modal contrastive pre-training for representing and traversing chemical space.  _ChemRXiv_ 10.26434/chemrxiv–2023–bdkgm (2023).
  * [79] VandeVondele, J. & Hutter, J.  Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases.  _J. Chem. Phys._ 127, 114105 (2007).
  * [80] Goedecker, S., Teter, M. & Hutter, J.  Separable dual-space Gaussian pseudopotentials.  _Phys. Rev. B_ 54, 1703–1710 (1996).
  * [81] Sun, J., Ruzsinszky, A. & Perdew, J. P.  Strongly constrained and appropriately normed semilocal density functional.  _Phys. Rev. Lett._ 115, 036402 (2015).
  * [82] Hutter, J.  https://github.com/juerghutter/GTH/ (2021).
  * [83] Miceli, G., Hutter, J. & Pasquarello, A.  Liquid water through density-functional molecular dynamics: Plane-wave vs atomic-orbital basis sets.  _J. Chem. Theory Comput._ 12, 3456–3462 (2016).
  * [84] Yao, Y. & Kanai, Y.  Free energy profile of NaCl in water: First-principles molecular dynamics with SCAN and $\omega$B97X-V exchange-correlation functionals.  _J. Chem. Theory Comput._ 14, 884–893 (2018).
  * [85] Martyna, G. J., Klein, M. L. & Tuckerman, M.  Nosé–Hoover chains: The canonical ensemble via continuous dynamics.  _J. Chem. Phys._ 97, 2635–2643 (1992).
  * [86] Guidon, M., Hutter, J. & VandeVondele, J.  Auxiliary Density Matrix Methods for Hartree-Fock Exchange Calculations.  _J. Chem. Theory Comput._ 6, 2348–2364 (2010).
  * [87] Pagotto, J., Zhang, J. & Duignan, T. T.  Predicting the properties of salt water using neural network potentials and continuum solvent theory.  _ChemRXiv_ 10.26434/chemrxiv–2022–jndlx (2022).
  * [88] Hockney, R. W.  _Computer Simulation Using Particles_ (A. Hilger, 1988).
  * [89] Plimpton, S.  Fast parallel algorithms for short-range molecular dynamics.  _J. Comput. Phys._ 117, 1–19 (1995).
  * [90] Spivak, M. _et al._ VMD as a Platform for Interactive Small Molecule Preparation and Visualization in Quantum and Classical Simulations.  _J. Chem. Inf. Model._ 63, 4664–4678 (2023).
  * [91] Hess, B., Holm, C. & van der Vegt, N.  Modeling Multibody Effects in Ionic Solutions with a Concentration Dependent Dielectric Permittivity.  _Phys. Rev. Lett._ 96, 147801 (2006).
  * [92] Wolfram Research Inc.  Mathematica.  Wolfram Research, Inc. (2019).
  * [93] Rasaiah, J. C. & Friedman, H. L.  Integral equation methods in the computation of equilibrium properties of ionic solutions.  _J. Chem. Phys._ 48, 2742–2752 (1968).
  * [94] Vrbka, L. _et al._ Ion-specific thermodynamics of multicomponent electrolytes: A hybrid HNC/MD approach.  _J. Chem. Phys._ 131, 154109 (2009).
  * [95] Herrington, T. M., Pethybrldge, A. D. & Roffey, M. G.  Densities of Aqueous Lithium, Sodium, and Potassium Hydroxides from 25 to 75 °C at 1 atm.  _J. Chem. Eng. Data_ 31, 31–34 (1986).
  * [96] Sipos, P. M., Hefter, G. & May, P. M.  Viscosities and densities of highly concentrated aqueous MOH solutions (M+ = Na+, K+, Li+, Cs+, (CH3)4N+) at 25.0°C.  _J. Chem. Eng. Data_ 45, 613–617 (2000).
  * [97] Lide, D. R. (ed.) _CRC Handbook of Chemistry and Physics, Internet Version_ (Taylor and Francis, Boca Raton, FL, 2007).
  * [98] Pitzer, K. S. & Mayorga, G.  Thermodynamics of Electrolytes. II. Activity and Osmotic Coefficients for Strong Electrolytes with One or Both Ions Univalent.  _J. Phys. Chem._ 77, 2300–2308 (1973).
  * [99] Kusalik, P. G. & Patey, G. N.  The thermodynamic properties of electrolyte solutions: Some formal results.  _J. Chem. Phys._ 86, 5110–5116 (1987).
  * [100] Yeh, I.-C.  System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions.  _J. Phys. Chem. B_ 108, 15873–15879 (2004).