Ab Initio Molecular Dynamics Simulations of the Influence of Lithium Bromide on the Structure of the Aqueous Solution-Air Interface.

We present the results of ab initio molecular dynamics simulations of the solution-air interface of aqueous lithium bromide (LiBr). We find that, in agreement with the experimental data and previous simulation results with empirical polarizable force field models, Br- anions prefer to accumulate just below the first molecular water layer near the interface, whereas Li+ cations remain deeply buried several molecular layers from the interface, even at very high concentration. The separation of ions has a profound effect on the average orientation of water molecules in the vicinity of the interface. We also find that the hydration number of Li+ cations in the center of the slab Nc,Li+-H2O ≈ 4.7 ± 0.3, regardless of the salt concentration. This estimate is consistent with the recent experimental neutron scattering data, confirming that results from nonpolarizable empirical models, which consistently predict tetrahedral coordination of Li+ to four solvent molecules, are incorrect. Consequently, disruption of the hydrogen bond network caused by Li+ may be overestimated in nonpolarizable empirical models. Overall, our results suggest that empirical models, in particular nonpolarizable models, may not capture all of the properties of the solution-air interface necessary to fully understand the interfacial chemistry.

Introduction

From seawater to the cellular environment and in many industrial processes, ionic salts are a major component of aqueous systems. Great progress has been made in understanding the influence of monatomic salt ions on, for example, the air–liquid interface. The determination of the tendency of larger anions such as Br– and I– to reside at or near the edge of an aqueous interface, in opposition to classical solution theories,[1] resulted from an intensive combination of experimental, simulation, and theoretical work over the course of many years.[2−12] The molecular simulation aspects of this story relied on the ability of classical polarizable models to describe the ion–water interactions in a way that could account for the changing electronic properties of the ions, leading to a polarization-driven change in the ionic solvation and the appearance of a free energy minimum for ions near the interface.[3−5,12−14] However, a debate continues regarding whether reparameterized nonpolarizable models may still be able to reproduce the correct interfacial behaviors.[15−17] Polarizability may also be important for the correct description of ionic interactions at other interfaces, for example, binding of Ca2+ to phospholipid bilayers[18] or interactions of Na+ with silica nanopores.[19]

Given the tremendous interest in and importance of better understanding the air–water interface of ionic solutions, initially, it seems remarkable that simulation methods which go beyond empirical polarizable models have not been used much to study this system. It is only recently that progress in computational technology has made ab initio molecular dynamics (AIMD) tractable for the study of the interfacial properties of simple neat water.[20−24] The introduction of salt ions adds new challenges, including ensuring that the same methods and basis sets used for pure water are adequate when ions are added and considerably longer simulation times to converge structural and dynamical properties of interest. Two recent studies have investigated the air–liquid interface of aqueous sodium chloride using Car–Parrinello molecular dynamics (CPMD).[25,26] Previously, quantum chemical investigations of ion interactions in interfacial aqueous systems were limited to rather small clusters[27−29] or single ions in a water slab,[14,29] which while relevant, are not sufficient to describe a more concentrated solution with an extended interface.

In this study, we examine the effect of addition of lithium bromide (LiBr) salt to water using density functional theory (DFT)-based simulations, in particular AIMD, and more specifically Born–Oppenheimer molecular dynamics (BOMD). In BOMD, the DFT functional is fully minimized at each timestep, in contrast to the more efficient but somewhat approximate CPMD methodology. On the other hand, BOMD allows the use of longer integration timesteps than CPMD. Our choice of LiBr is motivated by our previous investigations using BOMD to study the properties of formic acid and other molecules at a water–air interface,[30,31] which in turn was an attempt to use simulations to interpret experimental investigations of molecular collisions with liquid microjets, and the probability that the molecules are absorbed or reflected by the liquid surface.[32,33] These experiments must use very high concentrations (∼8 M) of lithium bromide or some other highly soluble salt to avoid water evaporation. It seems plausible that these high salt concentrations may also impact on the interactions between incident molecules and the liquid surface. Concentrated aqueous solutions of LiBr are also used as liquid desiccators in industrial applications,[34] and so, their interfacial properties may be of interest in these fields as well.

We have simulated aqueous LiBr in a pseudo-two-dimensional slab geometry, focusing on the interfacial structure of the solutions, but our results are also relevant to answer some basic questions about both Li+ and Br– solvation in bulk solution. In particular, resolving the Li+ hydration number has been a long-standing issue. All simulations with nonpolarizable empirical models predict a stable tetrahedral coordination shell with four hydrating water molecules, and no otherwise reliable nonpolarizable empirical model for the Li+–water interaction has been shown to predict much deviation from this tetrahedral hydration.[35] One interpretation of this very stable tetrahedral coordination is that the Li+ ion acts in some ways as a larger ion which includes the four hydrating water molecules.[17,36] Available simulation data using AIMD or hybrid QM-MM methods also indicate a tetrahedral coordination for Li+;[37−39] however, these studies did not include dispersion effects in their computations of the quantum mechanical potential energy surface, and because inclusion of dispersion has been shown to be critical to describe the structure of pure water accurately,[23,24,40−43] it is likely that dispersion will be important in ionic solutions as well. Inclusion of dispersion has already been found to be important in CPMD simulations of the solution–air interface of aqueous sodium chloride.[26] Recent quantum chemical results also demonstrate that the tetrahedral cluster of Li+·4H2O has higher binding energy per molecule than the hexagonal Li+·6H2O cluster but with some evidence for a more labile solvation shell than that predicted by empirical methods.[44]

By contrast, experimental data derived from neutron scattering in LiCl aqueous solutions measure significantly larger hydration numbers, about 4.3 at higher concentration and as high as 5.9 as the ion concentration is lowered.[45−47] Available simulation data using empirical polarizable models[48,49] indicate that they can also be parameterized to generate Li+ coordination numbers above four, which are in better agreement with the experimental data. The simulations we present in this work using a fully quantum mechanical description of the potential energy surface and fewer adjustable parameters may represent a significant step forward in more accurately describing the Li+ coordination shell.

Computational Details

For initial equilibration and for comparison with AIMD, we used LAMMPS[50] (11 Aug 2017 version) to run MD simulations using empirical nonpolarizable models. The TIP4P-Ew model[51] was applied for water atoms along with the associated Joung–Cheatham ion models[52] designed for use with this water model. In the two-dimensional (2D) periodic slab geometry, we applied the correction to the three-dimensional (3D) Ewald summation including the usual fictitious extra space 3L added to prevent interactions between periodic images in the z direction,[53] where the z direction is defined to be perpendicular to the interface. Walls with repulsive Lennard-Jones potentials were placed at z = 0 and z = L in order to prevent any rare evaporating molecules from being lost in the nonperiodic z dimension. A Langevin thermostat was applied to maintain the temperature T near 300 K.

After initial equilibration over several nanoseconds with the empirical models, molecular snapshots were used to instantiate ab initio BOMD simulations. These simulations were run using the QUICKSTEP module implemented in the CP2K code[54] with a timestep of 0.25 fs. The energy density cutoff used was 280 Ry, the default value in CP2K. For the density functional computations of the potential energy, the BLYP exchange–correlation functional[55,56] was used in conjunction with Grimme’s D2 dispersion correction.[57,58] Available results comparing the effect of Grimme’s more advanced D3 correction with D2 on aqueous systems indicate that the main impact is an improvement in modeling dynamical properties such as the diffusion constant.[59] Because dynamical properties are not the focus of this work, we used the D2 correction to allow direct comparison with our previous investigations of similar interfacial aqueous systems.

The DZVP-MOLOPT-SR-GTH basis functions were used for the basis sets along with the BLYP–GTH pseudopotentials.[60] This combination of the level of DFT and basis set has been shown to be sufficient to describe both bulk water and the air–water interface with a high level of accuracy;[23,42,59] however, the equilibrium water density at T = 300 K and zero pressure is somewhat overestimated.[30,41] There are results showing that a triple-zeta basis set (TZV2P) may better reproduce the correct density of water at T = 300 K;[24] however, this larger basis set would have significantly increased the computational cost of our simulations. Instead, we have run our simulations at T = 330 K in order to more closely match the expected density at room temperature. The temperature was held constant by using a Nose–Hoover chain thermostat of length 3 and time constant 50 fs applied on the massive degrees of freedom.[61,62] The BOMD simulations were run in a 3D periodic geometry with sufficient empty space to prevent significant errors because of interactions between periodic images in the z direction.

We found that the potential energy of the AIMD simulations converged after only a few ps of simulation time, after which, we began to accumulate data. A snapshot of the equilibrated AIMD simulation system with 44H2O molecules and 10 LiBr pairs is shown in Figure . We define the z axis consistently as the direction perpendicular to the liquid–air interface. In all of our data analyses, we set the origin of the z axis at the Gibbs dividing surface, defined as the point where the number density of the liquid has dropped to half of the density in the center of the slab, as shown schematically in Figure . Furthermore, we take advantage of the symmetry of the slab to average results on either side of the center of mass and show only one half of the slab in our figures.

Figure 1

Snapshot of the AIMD simulation with 44H2O molecules and 10 LiBr ion pairs (system SHC in Table ).

Table 1

Summary of All Simulation Runsa

label	nH2O	nLi+,Br–	mass % LiBr	Lx,y/Å	Lz/Å	trun	Nrun
Empirical Models
SW	64	0	0	12.5	60.0	20 ns	2
SLC	58	3	20.0	12.5	50.0	20 ns	4
SHC	44	10	52.3	12.5	50.0	20 ns	9
MW	125	0	0	16.5	50.0	20 ns	2
MLC	117	4	14.1	16.5	50.0	20 ns	6
MHC	85	20	53.1	16.5	50.0	20 ns	7
MHC-LZ	85	20	53.1	12.5	80.0	20 ns	4
LLC	480	16	13.8	23.0	80.0	10 ns	8
LHC	352	80	52.3	23.0	80.0	10 ns	6
DFT
SW	64	0	0	12.5	40.0	150 ps	1
SLC	58	3	20.0	12.5	40.0	150 ps	4
SHC	44	10	52.3	12.5	40.0	125 ps	4
MHC	85	20	53.1	16.5	40.0	150 ps	3

Column 1 is the label introduced for a reference in the rest of the article; the first letter indicates the system size, whereas the second and third letters indicate the salt concentration. Columns 2–4 list the number of water molecules, LiBr ion pairs, and the LiBr mass fraction, respectively. Columns 5 and 6 give the simulation box dimensions in Å, where 1 Å = 10–10 m. Column 7 gives the total time for each simulation, and column 8 is the number of independent simulations for each system.

Results and Discussion

We have run simulations with a variety of system sizes and geometries, as well as variations in the concentration of LiBr. In Table , we summarize all of the systems we have studied in this work, along with the simulation times for each. Labels are introduced for easy identification of each simulated system in the remainder of the paper.

The diffusion of the heavy Br– ion is comparatively slow, and therefore, long simulations were needed to ensure well-converged results. To somewhat mitigate the convergence problem, we have initiated simulations from several independent starting configurations in order to obtain reliable data as an average. Data from the AIMD simulations still have some significant error bars, which we show for some representative points in our plots. Block averaging of some of our data sets also demonstrated to our satisfaction that the simulations were well-equilibrated and converged (results not shown). The data for the empirical models are better converged, with relative errors generally less than 1%, and so, to aid the readability of the graphs, we do not show their error bars.

Size Dependence of Simulations in a Slab Geometry

Given the small size of the systems amenable to ab initio simulations, it is important to determine whether these systems are sufficiently large to accurately represent a real liquid–vapor interface. Available results with small systems (ca. 1–200 molecules) using AIMD to simulate pure water interfaces seem to be reliable.[23,24] When additional components are introduced, this question must be re-examined.

We have made use of simulations with empirical models to examine small systems with 64 or 125 molecules and compare the results with larger systems more commonly used to simulate interfacial systems with empirical methods. We show the number density histograms for water and each ionic species using empirical models in Figures and 3. There is some size dependence, in that the peaks due to molecular layering at the interface are somewhat exaggerated in the smaller systems, and there are still some density fluctuations in the middle of the slab. It is somewhat surprising that increasing the width of the slab in system MHC-LZ did not significantly change the height of the density peaks closest to the interface; rather, it seems that increasing the periodic box lengths L is more important. More prominent layering in systems with smaller x and y dimensions could be due to smaller fluctuations of the surface imposed due to the smaller system size. Overall, we find that the empirical force field results in the smallest systems we studied (SLC and SHC) are similar enough when compared with the larger systems, and assuming that the size dependence in AIMD simulations is not too different, we can view the DFT results in small systems as accurate.

Figure 2

Size dependence of the component density profile for simulations of concentrated LiBr solutions using empirical models. These (and all other similar histograms) were generated by averaging the number density over the course of the simulation within thin slices of the system determined by each atom’s distance z from the slab center of mass. Left: H2O (solid lines). Right: Ions only, Li+ (solid) and Br– (dashed).

Figure 3

Size dependence of the component density profile for simulations of dilute LiBr solutions using empirical models. Left: H2O (solid lines). Right: Ions only, Li+ (solid) and Br– (dashed). Note that the mass fraction of LiBr is different for systems MLC and LLC (∼14%) and for system SLC (20.0%).

Mass Density of Aqueous LiBr

The experimental mass density ρ of aqueous LiBr has been determined at a series of temperatures and compositions.[63] An empirical equation fit to these data gives estimates of ρ = 1.57 and 1.17 g/cm3 at LiBr mass percentages of 52.3 (system SHC) and 20.0 (system SLC), respectively, and T = 300 K. In Figure , we show the mass density profile from our simulations, as a function of the z position in the slab with respect to the interface.

Figure 4

Total mass density of aqueous LiBr solutions. Left: Empirical models Right: DFT-based AIMD. Error bars shown (in this figure and in others) are one standard error in the mean, derived from averages over several runs and both sides of the symmetric liquid slab.

The combination of empirical models we have used matches the experimental mass density very well, even with large ion concentrations. The match to the AIMD data is less clear, because of a combination of rather large error bars and the aforementioned molecular layering which complicates the analysis of small interfacial systems, but estimates of 1.7 ± 0.2 g/cm3 for systems SHC and MHC and 1.2 ± 0.1 g/cm3 for system SLC seem reasonable. While somewhat larger than the experimental measurement, within the statistical errors of our data, the results are in agreement. We note again that a mild overestimate (∼5–10%) of the liquid water density has been observed previously in AIMD simulations using methodology similar to that we have employed in this article.[30,41] The CPMD simulations using the BLYP-D2 method for a slab of concentrated aqueous NaCl run at T = 300 K also overestimate the mass density by a significant amount.[26] By running simulations at T = 330 K, we have somewhat corrected for this discrepancy, but there is still room for improvement.

At first glance, the large fluctuation in the mass density at the interface may seem surprising, compared with the much lower level of fluctuations seen in simulations of pure water liquid–vapor interfaces[22] or a low concentration aqueous sodium chloride interface.[26] Two main factors are worth considering here. First, the high concentration of salt naturally leads to large effects on the relative concentrations of different species. Second, as regards the mass density in particular, the large mass difference between Li+ and Br– causes the already significant fluctuations in the individual component number density profiles to be magnified in the total mass density profiles.

Comparison of Empirical Models with AIMD Results: Ion Density Profiles

In Figure , we show the average number density of all species as a function of z coordinate across the slab, computed in our DFT-based AIMD simulations, whereas in Figure , we compare the ion number density profile across the slab from both the empirical model simulations and the AIMD simulations. All of the simulations display some structuring in the ion density, which is most pronounced near the interface; the chief difference we see between the different models is a large degree of charge separation. With the empirical models, the peaks in both the Li+ and Br– concentrations are close to the same z position. In fact, the empirical models have Li+ slightly closer to the interface than Br–. Slight preference of Li+ for the interface has been noted by others in empirical model simulations[17] and interpreted as the tight tetrahedral coordination shell causing the Li+ to behave effectively as a much larger ion.

Figure 5

Density profiles of all species in AIMD simulations of LiBr solutions.

Figure 6

Density profiles of ionic species in highly concentrated LiBr solutions, comparing empirical models (black, green) with AIMD simulations (red, blue). Li+ number densities are the solid dashed lines. Br– densities are the dashed solid lines.

By contrast, in the AIMD simulations, the Br– ions tend to lie much closer to the interface than the Li+ ions. Our AIMD results are qualitatively similar to those that have been seen previously in simulations with polarizable models[4,12] and in experimental measurements[7,8] which show that larger anions such as Br– prefer the surface. In CPMD simulations of aqueous NaCl, only a small surface preference for Cl– versus Na+ is seen.[26]

Interfacial Orientation of Water Molecules

It is well known that most of the water molecules at an air–water interface, or indeed more generally at any hydrophobic interface,[64] reorient in such a way that the majority of the molecules have their dipoles parallel to the interface,[5,65] with a small bias into the bulk. In this way, the water molecules are able to maintain a maximum number of hydrogen bonds, just under three, rather than two as one might expect in the absence of reorientation.[66,67] Some alteration in the orientational tendencies of the water molecules might be expected because of the influence of the ions.

The average orientation of the water molecule dipole moments, where μH is defined as the vector from oxygen through the midpoint between two hydrogens, is shown in Figure . The average orientation of μH with respect to the z axis is affected by the z position as well as by the ionic strength. In particular, the charge separation in the AIMD simulations because of the different interfacial populations of the two ionic species, shown in the previous section (Figure ), causes an increased tendency for the water molecules in two molecular layers below the interface to have their dipoles projected somewhat out of the liquid slab. Empirical models also predict some changes in the orientational profiles because of the salt but not the large alteration in the dipole direction in the region just below the interface seen in the AIMD simulations.

Figure 7

Average orientation of water molecular dipole moments μH with respect to the z axis, as a function of location in the slab. Systems with 64 total molecules are on the left and with 125 total molecules on the right. Circles and solid lines are DFT results, and ×’s and dashed lines are empirical model results. Lines are cubic spline fits between the data points.

Direct comparison with the CPMD results for aqueous NaCl is difficult because only data for the average water orientation in the entire interfacial region (roughly |z| < 1.0 Å) are available,[26] showing only a small bias (∼−0.1) for water dipoles to orient into the slab. Our AIMD results show that even in the interfacial region proper, a small bias (∼0.1) remains for water dipoles to point out of the slab. It is likely that the large charge separation between Li+ and Br– is responsible for this difference.

Effects of Ions and the Interface on Inter-Molecular Bonding

Several different intermolecular noncovalent bonds are seen in our system. In addition to hydrogen bonding between water molecules, strong attractive interactions exist between Li+ cations and water oxygens and between Br– anions and water hydrogens, as well as ion-pairing interactions between Li+ and Br–. We use simple geometric criteria to define the existence of intermolecular bonds,[68] which have been validated and used in many previous studies. Hydrogen bonds between water molecules are defined as existing if the interoxygen distance rOO is less than the first minimum in the radial distribution function (RDF), that is, rOO < 3.5 Å at lower ion concentrations and rOO > 3.7 Å at higher concentrations, at the same time as there is an intermolecular oxygen–hydrogen distance rOH < 2.5 Å and the O···H–O angle <30°. Ion–water bonds are also defined based on the first minimum in the RDF, rOLi < 2.5 Å and rHBr < 3.1 Å. Finally, an ion pair exists if rLi < 3.3 Å.

Coordination Numbers of Li+ and Br–

We have plotted the average ion coordination numbers as a function of z, ⟨Nc⟩(z) in Figures and 9 for Li+ and Br–, respectively. Hydration numbers for ion–water contacts are shown, in addition to the average number of ion–ion pairs and the total coordination numbers, including both.

Figure 8

Average number of coordinating molecules per Li+ ion, ⟨Nc,Li⟩(z) in both DFT (left) and empirical model (right) simulations. Coordination to water molecules (hydration number) is shown in black, while coordination to Br– and the total coordination numbers are shown in blue and red, respectively.

Figure 9

Average number of coordinating molecules per Br– ion, ⟨Nc,Br⟩(z) in both DFT (left) and empirical model (right) simulations. Coordination to water molecules (hydration number) is shown in black, while coordination to Li– and the total coordination numbers are shown in red and blue, respectively.

The Li+ hydration numbers as predicted by the AIMD simulations vary somewhat according to the position of the ions in the slab, but in the center of the slab at low concentration, it is about 4.8. At low ion concentration, the degree of ion pairing is negligible, while at higher concentration, most of the ions are paired with a single counterion. The total coordination number at high concentration is approximately 5.2, with the hydration number alone somewhat lower because of the ion pairing. The coordination number becomes lower for the rare Li+ ions near the interface.

The empirical model simulations, however, show that the Li+ ion invariably coordinates to a total of four other molecules, with the identity of these molecules varying according to the ion concentration. It is clear that the lack of variation in the total Li+ coordination number is somewhat unusual, suggesting an overly structured tetrahedral coordination shell with little exchange of solvent molecules with the bulk.

The degree of agreement between simulations and the neutron scattering experiments regarding the Li+ coordination shell has been subjected to some debate. Recent estimates of the hydration number from neutron scattering range between 4.8 and 5.9 at low concentration down to approximately 4.3 at higher concentration.[46,47] While our simulation results deviate from the most recent data showing a hydration number close to 6,[47] it is significant that we find hydration numbers well above the result of 4.0 seen in all nonpolarizable empirical models.[35]

As for the Br– hydration number, there is considerably closer agreement between our two methodologies. At high concentration, in the center of the slab, the hydration number is about 6.0 in the empirical model simulations, rising slightly to about 6.5 in the DFT-based simulations. Both methods predict some ion pairing at high concentration. The AIMD results show some variation in the ion pairing, as the average number of ion pairs drops from around 1 in the center to below 0.5 near the interface. This is consistent with how Li+ tends to avoid the interface in the DFT simulations. The hydration number at low concentration, by contrast, is somewhat lower in the AIMD case (∼6.0) compared with the empirical models (∼6.7). All of the simulations show a gradual decrease in the Br– coordination numbers as the ion is closer to the interface.

Hydrogen Bonding and Water–Ion Contacts

In Figure , we plot the average number of intermolecular bonds formed between water and other species. Overall, the two different simulation models give similar results. In agreement with many previous results,[66,67] pure water maintains just under four hydrogen bonds in the bulk of the slab, gradually lowering by about one H bond to just under three at the interface, as water molecules reorient to maximize the number of bonds they can form.

Figure 10

Average number of intermolecular bonds per water molecule ⟨nHB⟩(z) in both DFT (left) and empirical model (right) simulations. Water–water hydrogen bonds are shown separately as solid lines. The total number of intermolecular bonds including water–ion bonds is shown as dashed lines.

As the ion concentration increases, some of the water–water hydrogen bonds are replaced by water–ion contacts, but the total number of molecules in each water molecule’s coordination shell remains rather constant, especially in the AIMD simulations. In the empirical model simulations, we see a significant reduction in the average number of coordinating molecules as concentration increases. At the interface (z = 0), the total number of coordinating molecules per each water molecule drops from ∼2.7 in pure water down to ∼2.2 in highly concentrated LiBr. A plausible explanation for this large effect is that the overly stable tetrahedral Li+ coordination shell, combined with the closer approach of Li+ to the interface, causes a significant disruption to the overall hydrogen bond network in general. By contrast, it would appear that the more labile Br– solvation shell, which must influence the interface more than Li+ in the AIMD simulations, does not perturb the hydrogen bond network as drastically, and even the rare Li+ cations which may approach the surface do not play as large a role as they do in the empirical simulations.

The CPMD simulations of the aqueous NaCl interface showed ⟨nHB⟩ ≈ 2.8 water–water hydrogen bonds to be maintained in the center of the slab, lowering to ∼2.0 in the interfacial region.[26] These results, for an ion concentration (5.3 M) in between our low and high concentration systems, are in fair agreement with ours, which show ⟨nHB⟩ ≈ 2.8 and ∼2.3 for bulk and interfacial water–water bonds, respectively, at low concentration, and ⟨nHB⟩ ≈ 1.5 and ∼1.1 at high concentration. As one might expect, the larger Br– anion participates in more intermolecular water–ion bonds than Cl–, leading to larger reductions in the total number of water–water bonds.

Conclusions

We have completed a comprehensive AIMD study of a system of aqueous lithium bromide in a pseudo-2D geometry, allowing us to investigate the interfacial properties of the solution. Despite the size limitations inherent in ab initio simulations, we are confident that our methodology provides an accurate representation of the real air–liquid interface of aqueous LiBr.

We find that, in agreement with the experimental results and empirical model simulations with polarizable force fields, Br– anions have a much higher affinity for the interface compared with Li+ cations. This large separation between the ionic charges in molecular layers just below the interface leads to water molecules affected by the charge separation now preferring to orient their dipoles toward the surface, rather than slightly into the bulk as in pure water or in empirical simulations regardless of the salt content. It remains to be seen if this change may affect the physicochemical properties of the interface such as the adsorption and absorption propensities of different incident molecules.

Our results predict Li+ hydration and/or coordination numbers significantly above 4. Nonpolarizable empirical models have had difficulty predicting the correct structure of the Li+ solvation shell, tending to predict overly stable tetrahedral coordination. Our AIMD results are in better agreement with the most recent experimental data for the Li+ hydration number. The combination of better modeling of the Li+ solvation shell and more accurate prediction of the ion density profile suggests that nonpolarizable empirical model simulations might overestimate the reduction in the number of hydrogen bonds maintained by water molecules near the interface.

We have compared AIMD results with nonpolarizable empirical models. A worthwhile continuation of this study would be to compare the results of empirical polarizable models based on Drude oscillators[69] or other methodologies to investigate exactly how good the qualitative agreement we have noted is, in particular as regards Br– ion’s preference for the surface. It would also be interesting to examine the Li+ solvation shell with other analyses, including, for example, computing potentials of mean force with respect to the interface position, or a more rigorous analysis of the coordination shell geometries to understand the balance between tetrahedral and hexagonal solvation structures or other structures with different geometries. This would be a viable direction to pursue in future work.

One of the interesting applications for this work will be in the area of reactions in and on atmospheric aerosols, especially in a marine environment where salt ions are a major component. Aqueous sodium chloride was studied previously using computational methods similar to what we have used in the current paper.[25,26] Comparison between our study using LiBr and the previous work on NaCl demonstrates that our methodology should be equally applicable to aqueous interfaces including all kinds of alkali halide salts. Including salt in studies of reactions on aqueous surfaces could be an important part of understanding important interactions of gas-phase molecules with atmospheric and marine aqueous interfaces.

The focus of this work has been on structural quantities. Dynamical properties such as diffusion or other transport properties, or rotational correlation functions, are also of interest, and have been much studied in similar interfacial systems. However, the computation of transport properties, or indeed any time-dependent properties, raises challenging questions such as how to decide if a given ion or molecule should contribute to the average of a quantity computed in a particular region of the interface as it moves over the course of the simulation.[70] In addition, some quantities such as the diffusion constant may be anisotropic in an interfacial geometry.[71] These would be interesting questions to explore further in subsequent work.

Our results show that BOMD simulations based on a fully ab initio DFT potential energy surface can accurately describe many of the challenging bulk and interfacial properties of aqueous alkali halide solutions. Our future work will build on this study by introducing interactions between the solution and other molecules of interest, both at the interface and in the bulk.