The structure of concentrated Li-ammonia solutions as derived from MD simulations ∗

The results of Molecular Dynamics simulations of lithium-ammonia solutions over the whole concentration range from 0.5 to 19.6 MPM at 240 K are reported. The pseudopotential theory is employed at the higher concentrations and the direct contribution to the total potential which has been derived from ab initio calculations has been supplemented by three-body terms. The resulting partial structure and radial distribution functions are compared in detail with recent X-ray and neutron diffraction studies with isotopic substitution. There is an overall good agreement between simulation and experiment. Differences are discussed. The solvation number of the lithium ion is found to be six for the dilute solution and four for the higher concentrations in analogy to the hydration of Li found for various aqueous solutions.


Introduction
Recent investigations of various metal-ammonia solutions by Molecular Dynamics (MD) simulations demonstrated the importance of three-body contributions for the structure of the solvation shells of the ions [1][2][3][4].In addition, the comparison with neutron and X-ray diffraction studies showed some discrepancies, especially with the latest neutron diffraction studies with isotopic substitution (NDIS) at high concentrations [5][6][7][8][9].Although earlier simulations of very dilute Li-ammonia solutions did not show any effect of three-body contributions to the lithium ion solvation [10], we decided to simulate again Li-ammonia solutions over the whole concentration range where the pair potentials had been supplemented by three-body terms.
By applying the pseudopotential theory, effective renormalized interatomic potentials were derived through which electrons can be excluded from an explicit consideration and classical computer simulation methods can be used to calculate the structural and dynamical properties of concentrated metal-ammonia solutions [11].Thus, the total potential which describes the effective interactions between the sites of i and j in the metal-ammonia system containing free unlocalized electrons is given by where r is the distance between the two sites belonging to different molecules, between two metal ions, or between a metal ion and an ammonia molecule site.The direct contribution is based on ab initio calculations while the indirect contribution can be derived solely from the pseudopotential theory provided that the metal concentration is sufficiently high for the electrons not to be localized.In this way, the structure of a 11.64 MPM Li-ammonia solution (about half saturated) has been derived from an MC simulation with a rigid ammonia model [12].The reliability of the thus achieved results could not be checked because of a lack of detailed experimental information for this concentration.
In this work, we report on the simulations of Li-ammonia solutions over the whole concentration range.In the next chapter we discuss the potentials and the details of the simulations followed by a chapter where the structure of such solutions and its dependence on the lithium concentration is discussed in comparison with the results from neutron and X-ray diffraction studies.

Potentials and details of the simulations
A flexible four site model was employed to describe the ammonia-ammonia interactions.The experimental gas phase geometry, with an N-H distance of 1.0124 Å and an HNH angle of 106.67 • , was taken from reference [13].The partial charges of 0.2674|e| on the hydrogen and −0.8022|e| on the nitrogen atom are taken from SCF calculations of the ammonia molecule [14].
All of the direct contributions to the potentials have been derived from ab initio calculations and are given for the ammonia-ammonia interactions by [15]: and for the Li-Li and Li-ammonia interactions by [12]: where the energies are given in kJ/mol with r in Å.For further details of the potentials and especially the indirect contributions [11], the reader is referred to the references cited.
The results of an MD simulation of a saturated Li-ammonia solution with these potentials has led to some discrepancies with the experimental findings.As a consequence of these discrepancies we have modified the potentials by introducing a 3-body term for the Li-ammonia interactions as described in detail in reference [10], although the simulation of a single lithium ion in 201 ammonia molecules did not significantly change the Li-N radial distribution function (RDF).The 3-body term employed is given by [10]: where r 1 and r 2 are the distances between Li + and each of the two nitrogen atoms and r 3 is the distance between the two nitrogen atoms.The fitting parameters are found to be: a = 656.92kJ/mol, b = 0.2261 Å−1 , and c = 0.9645 Å−1 .There was no need to introduce a 3-body term for the ammonia-ammonia interactions as discussed in reference [16].With these potentials four different MD simulations have been performed, the details of which are collected in table 1. Due to the strong screening of the Coulombic interactions by the indirect contributions to the total potential (see reference [12]) there was no need to use the Ewald summation.Instead, the shifted force method with a cut-off length of half the box size has been employed together with periodic boundary conditions.The time step length was chosen to be 0.125 fs.

Radial distribution and structure functions
The site-site radial distribution functions (RDFs) together with the running integration numbers are shown in figure 1 for four different lithium ion concentrations according to table 1.This report presents exclusively the results for the simulations including 3-body interactions.The characteristic values of the RDFs are presented in table 2.   Peak splitting (see figure 1).
It can be seen from the Li-N RDFs that the first peak is most pronounced for the single ion simulation where not only the first but also the second solvation shell can be formed undisturbed by other lithium ions.The results for this dilute case reproduce the findings presented already in reference [15] and proved to be unchanged when 3-body interactions were included [10].In accordance with these findings the solvation number remained six while for the higher concentrations the addition of the 3-body term reduced the solvation number to four.This change leads to an agreement with the NDIS measurements of a saturated Li-ammonia solution [6].An experimental confirmation of a solvation number of six for the dilute solution does not exist.The hydration number of the lithium ion is discussed controversially as well.There are, similarly, indications for 6 at low concentrations and for 4 at higher ones (see e.g.[17]).
It should be mentioned here that the NDIS measurements show a prepeak at 1 Å−1 which is only indicated in the simulation data at medium concentrations and reduces to a shoulder at the saturated solution.This prepeak is expected to reflect a medium range structure in real space but such structures can neither be seen in the RDFs from the experiment nor in those from the simulation.There exists only a well defined second solvation shell in the very dilute solution (see figure 1 and table 2).There is -in the limits of uncertainty -no difference between simulation and experiment in the positions of the first peaks in the Li-N and Li-H RDFs.It should be noted that the first maximum in the Li-N RDF is higher for the saturated solution than that for the 6.4 and 11.8 MPM ones, indicating a more pronounced solvation shell, which might indicate the cluster formation with the increasing concentration.
Similar conclusions can be drawn from the Li-H RDFs.Here, as well as in the Li-N case the first neighbor distances are slightly larger for the dilute versus concentrated solutions (see table 2).While for the 6.4 and 11.8 MPM solutions both the Li-N and Li-H RDFs very slowly but monotonously approach the uniform distribution, there is an indication of a slightly further ranging structure for the saturated solution.
With the screening of the Coulomb interactions in the pseudopotential theory, obviously the 3-body term in the Li-N and Li-H interactions becomes significantly more important than in the single ion case where the electron is simply neglected.This is very different from the simulation of the beryllium ion in water, where the 3-body term immediately reduced the hydration number from six to four even in the single ion case [18,19].
The relatively small number of Li + leads to a strong noise in the Li-Li RDF but this limitation does not hide the obvious effect that the ion-ion structure is significantly more pronounced for the 6.4 MPM solution than for the other two concentrations.This special feature is also reflected in all three solvent RDFs (see below).Consequently, in spite of the much smaller average ion density, the number of the nearest neighbors is about 0.7 while it is only about 0.5 for the 11.8 and 19.6 MPM solutions.The first neighbor Li-Li distance is about the same for all concentrations.The meaning of the Li-Li RDF for the solution structure is not yet clear.It deserves a further analysis to be derived from the simulation data.
The most detailed experimental confirmation on the change of the solvent structure with the lithium ion concentration has been derived from the neutron diffraction study with hydrogen isotope substitution [9].The measured partial structure functions (pstfs) and the RDFs for N-N, N-H, and H-H will be compared here with those derived from simulations with lithium ion concentrations of 6.4, 11.8, and 19.6 MPM.They are presented in figure 2. The comparison between simulation and experiment shows that the main peaks in all three pstfs and all concentrations are in the range 1.9-2.2Å−1 where for the highest concentration the experimental ones are slightly shifted by about 0.1 Å−1 to larger k.There is an agreement between simulation and experiment as far as the changes of position and height of the main peak with concentration are concerned.The positions shift in all cases to lower k-values with the increasing lithium concentration while its height decreases for the N-N pstf and increases for N-H and H-H with concentration.Simulation and experiment also show agreement in respect to the shoulder at about 2.5 Å−1 in the N-N pstf for the saturated solution.In addition, the second and the third small peaks at about 4 and 6 Å appear at similar positions.The pstfs N-H and H-H from the simulations practically disappear for k-values larger than 6 Å−1 , slightly different from the experimental results.The only significant discrepancy between simulation and experiment concerns the prepeak at about 1 Å−1 .It is most obvious and well pronounced in the total structure function for the saturated solution [9].In the simulation data, all three pstfs and at all concentrations also show a prepeak, but much less pronounced, which additionally reduces to a shoulder when the concentration is increased.It remains unclear why this discrepancy has been found.
According to the overall agreement even for the pstfs, a similar agreement between simulation and experiment for the partial RDFs is expected.The broad peaks resulting from the NDIS measurements [9] make it difficult to determine the positions and heights of the peaks in the RDFs and thereby the running integration numbers.For a better analysis, Gaussian functions have been fitted [8,9].The problems connected with this kind of analysis become obvious when the positions and the areas under the various peaks in the N-N RDF from the X-ray measurement [8] and the NDIS experiment [9] are compared.The same information is easily deduced from the simulation.The positions of the peaks and the number of the nearest neighbors (running integration up to the first minimum) for all partial RDFs and the different concentrations are given in table 2.
The position of the main peak in the N-N RDF is shifted to shorter distances by about 0.15 Å and its height decreases in going from the dilute to the saturated solution.The integration of the first peak up to about 5 Å (the first minimum) leads to a change in the number of the nearest neighbors from about 12 to 7.5 (figure 1 and table 2).This means a very good agreement with the experimental results, considering the difficulties in the analysis of the measured pstfs mentioned in the preceding paragraph [8,9].The N-N RDFs show a remarkable feature.The first and the second neighbor shells are most pronounced for the 6.4 MPM solution.It indicates -together with similar features in the N-H and H-H RDFs -a special solvent structure in the range of medium concentration which is not reflected in the experimental results and has yet to be analyzed more in detail.The medium distance range order in real space expected from the prepeak at about 1 Å−1 in the experimental N-N pstf is not reflected in the corresponding RDF [9].It might follow again from the difficulties in the analysis of experimental pstfs mentioned above.It remains unclear whether the shallow maxima and minima beyond 6 Å in the simulated N-N RDF lead to the weak prepeak found also in the simulated N-N pstf at 1 Å−1 (figure 2).
The characteristic data for the N-H RDF show -again in the limits of uncertainty -good agreement between simulation and experiment [9] for the dilute and the saturated solution.This is true not only for the position of the main peak but also for the number of the nearest neighbors.In both investigations a pronounced shoulder at about 2.5 Å is obvious which reduces with the increasing concentration.The formation of the shoulder has been discussed in reference [15] for the dilute solution.It is a result of the existence of two kinds of hydrogen atoms.The ones at the shorter distances are those which form bonds with the neighboring molecule.This is very similar to the water case where the integral over the first part of the double peak has often been employed as a measure of the number of hydrogen bonds formed.
But it has to be kept in mind that a distance criterion alone cannot be a correct measure of hydrogen bonding.The relative orientation of the two molecules has to be included in the definition.But having employed this estimation here would mean that the largest number of hydrogen bonds between solvent molecules is formed in the 6.4 MPM solution (see table 2 and reference 20) in accordance with the result for the N-N RDF mentioned above.The positions of the main peaks of the H-H RDFs again agree between simulation and experiment but the change of the RDF with concentration is much more detailed in the simulation when compared with the experiment.This is again especially true for the solution with the concentration of 6.4 MPM.

Structure of the solvation shells
The geometrical arrangement of the ammonia molecules in the solvation shells of the lithium ions can be deduced from the simulation by the calculation of the distribution of cos ϑ, where ϑ is defined as the N-Li-N angle.The result is depicted in figure 3, where only the ammonia molecules in the first solvation shell are included in the distribution.It is obvious from figure 3 that the six ammonia in the dilute solution are arranged octahedrally around the lithium ions.Integration over the distribution shows that there is another solvation shell molecule at about 180 • and four more at about 90 • .For all three higher concentrations the solvation number is four (table 2).The distributions of cos ϑ for these concentrations are very similar.There exists only one broad maximum around 102 • indicating not very well pronounced nearly tetrahedral arrangement of the ammonia molecules in the first solvation shells.
The orientation of the ammonia molecules in the solvation shells of the lithium ions can be described by the distribution of cos θ, where θ is defined as the angle between the dipole moment vector of the ammonia molecule and the vector pointing from the nitrogen atom towards the ion.The normalized distributions are depicted in figure 4.
For all four concentrations there is a preference for the dipole moment of the ammonia molecule pointing away from the ion.The distributions are quite similar for all three concentrated solutions but much broader than for the single ion case where solely the interaction between the solvation shell molecule and the ion determines the orientation because of the relatively weak interactions between the solvation shell molecules and those in the bulk.In the concentrated solutions, the solvation shells of different ions partly overlap and a competition between the two ions for the energetically most favorable orientation leads to a broad distribution.

Conclusions
The most important change caused by the introduction of three-body contributions is the reduction of the solvation number of Li + from 6 to 4 for the concentrated solutions.It remained 6 for the dilute one.This change of the solvation number with concentration in the ammonia solution is similar to what has been found for this ion in various aqueous solutions (see e.g.[17]).Very recent neutron diffraction studies with hydrogen isotope substitution [9] provide the possibility for a detailed comparison between simulation and experiment as far as the change of the solvent structure with the lithium ion concentration is concerned.There is an overall good agreement for the N-N, N-H, and H-H partial structure functions except for some differences in the prepeak near 1 Å−1 .But this discrepancy is not reflected in the partial radial distribution functions.
The N-N, N-H, and H-H RDFs from the simulations indicate the enhancement of the solvent structure in the medium concentration range, especially for the 6.4 MPM solution.This surprising feature is not obvious in the partial structure functions.Considering the overall good agreement between simulation and experiment, it is not to be expected that this special solvent structure is an artifact resulting from the simulation.Therefore, a further detailed analysis of the data and, most probably, additional simulations seem to be appropriate to clarify this interesting feature.

Figure 2 .
Figure 2. Partial structure functions for the N-N, N-H, and H-H interactions derived from MD simulations of lithiumammonia solutions for three different concentrations: 6.4 MPM (dotted), 11.8 MPM (dashed), and 19.6 MPM (solid).

Figure 3 .
Figure 3. Distribution of cos ϑwhere theta is defined as the nitrogenlithium-nitrogen angle (see insertion)calculated for the ammonia molecules in the first solvation shell of the lithium ion for four different concentrations: 0.46 MPM (dotted), 6.4 MPM (dash-dotted), 11.8 MPM (dashed), and 19.6 MPM (full).

Figure 4 .
Figure 4. Normalized distributions of cos θ for the ammonia molecules in the first solvation shell of the lithium ion for four different concentration: 0.46 MPM (dotted), 6.4 MPM (dash-dotted), 11.8 MPM (dashed), and 19.6 MPM (full).θ is defined in the insertion.

Table 1 .
Details of the simulations.

Table 2 .
Characteristic values of the radial distribution functions g ij (r) for Liammonia solutions at four different concentrations with the Li mole fractions c i : 0.0046, 0.064, 0.118, and 0.1958 denoted by 1, 2, 3, and 4, respectively.R l , R Ml , and R ml are the distances where for the lth time g ij (r) is unity, has a maximum, or a minimum, respectively.