1. Introduction

The large presence of the electrolyte solutions in our environment their inevitable utility in various biological, technical and industrial processes; and their preponderant presence in our life require us to take them into consideration and to study carefully their structures and their physicochemical properties, to understand their behaviors. They are, frequently, subjects of various areas of experimental and theoretical research; physics, chemistry, and biomedical.

For understanding the physicochemical properties and the reactions taking place in the electrolyte aqueous solutions, it is very important to conduct various structural studies of these solutions at the atomic level and in several conditions. In recent decades, much attention and great interest have been accorded to the interactions between water molecules and solvated ions to understand and explain the formation of solvation shells and the role of the hydrogen bonding. The aqueous solutions of LiCl present interesting properties which are studied by different methods at different concentration and thermodynamical states, so far many macroscopic properties of the aqueous LiCl solutions have been investigated [^{1 15}]. The main property is to form a glass through a metastable supercooled state when the temperature decreases [³,⁶,¹⁶,²⁰]. In this work, we investigate the hydration shell of the chlorine and lithium ions using a combination of the Hybrid Reverse Monte Carlo (HRMC) simulation with experimental scattering neutron data. The studied system is the LiCl6H₂O at the supercooled-liquid state (162 k). The obtained results by the Reverse Monte Carlo (RMC) in the latest paper [³] are improved and refined using an energetic constraint in the HRMC simulation, and they are used to highlight the hydrated structures of Li⁺ and Cl^-.

Although the relative strengths of the interactions between different ions and water molecules of aqueous electrolytic solutions in the first hydration shell are well investigated, the effect of ions on the H-bond between the first hydration shells is still discussed [⁵,¹³,¹⁶,³⁸]. In several papers, two situations have been proposed and used as a basic assumption to study the small and large ion hydrations; in one hand, the small ions interact strongly with the water molecules and break the H-bond, while the large ions interact weakly with water and facilitate the formation of H-bonds [²¹,²²]. Details of the simulation performed here are described in Sec. 2. In the following section, we present the obtained results and their discussions, and finally, in Sec. 4, we give a conclusion.

2. Computational Methods

The conventional algorithm of RMC uses the Markov process [²⁴-²⁷], but instead of minimizing the potential term as in the classical methods; molecular dynamics (MD) [²⁸] and Metropolis Monte Carlo (MMC) [²⁹] the difference between the calculated and the experimental data is the quantity to be minimized, Χ ², which is given by:

where G ^RMC(r _i ) and G ^EXP(r _i ) are the partial distribution functions obtained from the RMC configurations and experiment, and σ(r _i ) is an estimate of the experimental error. If Χ ² _new of a newly generated configuration is less than Χ _old; of the latest accepted configuration; the agreement between experimental and the current configuration is improved, and if Χ ² _new is increased, the new configuration is not rejected outright but accepted with the probability of acceptance, expressed by:

When a satisfactory agreement between experimental and theoretical data sets is obtained, detailed structural information such as the number and the average positions of the coordination and also the bond angle distribution functions can be calculated from the atomic networks. The idea is to obtain three-dimensional configurations in agreement with the available experimental data within a normal distribution of statistical errors.

The RMC method produces good fits to diffraction data, but with vastly different structures, and sometimes for ionic systems, the unrealistic structure appears in the pair distribution functions. McGreevy et al. [³⁰] pointed out that, in general, the lack of uniqueness is due to a deficiency in the three-dimensional information contained in the one-dimensional structure factor. Thus diffraction data constraints alone are insufficient for modeling physically realistic structures and possessing a thermodynamic aspect. To remedy these limitations, the most elaborate solution consists of the hybridization of RMC, and various types of constraints have been proposed to be applied in addition to the commonly geometrical constraints derived from the experimental data. The use of these constraints requires assumptions physically founded about the structure being modeled. Pikunic et al., [³¹] have used a bond angle constraint to model low density carbon, and O’Malley et al., [³²] have assumed a fully sp² bonded network to model glassy carbon.

In this work, we use a hybridization of the RMC algorithm using an energy constraint [⁴,⁶-³³-³⁶]. The combination is archived by adding a Boltzmann’s weighted energy term ΔU/K _B T to the total Χ ², in the same manner as other constraints are introduced; this is a hybrid of the Metropolis Monte Carlo (MMC) [²⁹] method and the RMC method. The agreement factor Χ ² becomes:

here, U denotes the potential energy penalty term, and ω is a weighting parameter and represents the temperature of the studied system. Acceptance criteria expressed by the conditional probability is now given as:

where ΔU = U _new - U _old is the change in energy due to the random Monte Carlo move, (Χ ² _new, U _new) and (Χ ² _old, U _old) are the agreement factors and energies of the new and old configurations, respectively. The energy of the system is calculated by using the Coulomb Buckingham potential expressed by:

where q _i , q _j are the charges of the individual ions, ε₀ is the permittivity of free space, σ _ij is a parameter characterizing the depth of the potential well, and ρ _ij is the slope of the short range exponential repulsion known as the hardness parameter [³⁷].

The Buckingham potential parameters are selected by fixing the hardness parameters for all the ions [³⁷]. This potential allows reducing the occurrence of very unlikely, high energy configurations, and consequently, limits the spectrum of acceptable configurations and eliminates unrealistic features. The atomic configuration used here is a cubic box of 31.0993 Å, corresponding to a number density of 0.09575 atoms/ų. It is composed of 2880 atoms: 864 Oxygens, 1728 Hydrogens, 144 Chlores, and 144 Lithium. The periodic boundary conditions are used to simulate the macroscopic level. Atoms of Oxygen (with a random orientation of the water molecule) form initially a face-centered cubic lattice, and the chlorine and lithium atoms are placed in the interstices.

3 Results and discussions

3.1 Partial correlation functions H(r)

In this work, The RMC modeling is conducted based on four experimental partial distribution functions (PDF’s): H ^EXPX _XH (r), H ^EXP _HH (r), H ^EXP HH(r), and H ^EXP _Clα(r), obtained from neutrons scattering experiment combined with the isotopic substitution [3]. Four types of correlations are described, where the subscript X defines all-atom species except the hydrogen one, while Cla represents the correlations between Cl and all the other species.

Figure 1 shows the experimental Partial Correlation Functions (PCF) H _ij (r) equivalent to the Partial Distribution Functions (PDF) G(r) (H(r) = G(r) - 1), and also those computed by RMC and HRMC for the supercooled-liquid state. The differences ΔH between them (H ^RMC _ij(r) - H _ij(r)) and (H ^HRMC _ij(r) - H ^EXP _ij(r)) are also drawn to permit an easy comparison.

Figure 1 Experimental and computed Partial Correlation Functions H _ij (r) for the supercooled-liquid state, ΔH _ij is the difference between them. 

All of the obtained results show a good agreement and a clear concordance with the experimental results. There is no discrepancy between RMC with or without the energy constraint; however, no conflict can be reported between the studied system and the used potential model. Consequently, the generated three-dimensional configurations are compatible with the experimental data, and thus, allow displaying all of the pair distribution functions. But it is noteworthy that the HRMC fit is better, and the corresponding ΔH is on average less than that of the conventional RMC.

3.2. Pair distribution functions g(r)

In this part, we observe the intern structure of water in the pure state at room temperature and in the solution through the water-water correlations, and then, we provide information on the hydration shells structure represented by the water-ion correlations, and finally, we discuss the ion-ion correlations. All of the water-ions and water-water pair distribution functions within the water structure and around ions (H/O - Li⁺, H/O - Cl, O - O, O - H, H - H ) are computed and to give supplement details to understand the structure of our solution, the coordination number for solute and solvent are also computed by using these functions. The running coordination number is obtained by integrating the respective pdfs, g(r).

n(r)=4πρ∫r2g(r)dr

where ρ is the number density, and r is the separation between the two species. The RMC simulation allowed us to generate three-dimensional configurations compatible with the experimental data, but we can notice that they present some anomalies in the form of artifacts or parasitic peaks accompanying the main peaks of coordinations. This can be due to the limited set of experimental data to only four PDFs and/or to small inconsistencies in experiment and/or to the non uniqueness problem [⁴]. The use of the inter-atomic energy term in HRMC is a very efficient tool against the local parasitic structures [⁴,⁶,³³,³⁵]. The obtained findings show a significant improvement of the pair distribution functions affected by artifacts, and we can observe the quasi-disappearance of the most important artifact peak located at 3.1 Å in the g _OO(r) (this artifact is marked by an arrow in Fig. 2a). Many other artifacts disappear, and the corresponding different pair correlation functions are suitably smoothed.

Figure 2 Pair distribution functions water-water: g _OO(r) (a), g _OH(r) (b), and g _HH(r) (c) at the supercooled-liquid state contrasted to the pure water at room temperature. 

3.2.1 Water structure

In the following passages, pair distribution functions that characterize the water-water correlations are discussed. To organize the discussions of our results and to make the analysis easier, the classification of the main peaks of the pair distribution functions shall be according to their natures; from the neighborhood of the water molecule; (the intermolecular interactions or hydrogen-bonding: H-bond) to its internal structure (the intramolecular correlations). The inter and intramolecular correlations of the water molecule are observed in Fig. 2, where all of the pair correlation of the molecule in the solution are contrasted to those of the pure water.

In pure water, molecules are usually arranged in a tetrahedral network of H-bonds, and this is manifested by the position of the first two peaks of the g _OO(r). The first one appears at r = 2.92 Å in the solution and at r = 2.8 Å in the pure water and a less intense second peak is spotted at the average distance of r = 4.5 Å in the pure water and at r = 4.7 Å in the solution (Fig. 2a). It is known that the second peak of the solution is affected by the interaction between ions and water molecules in the ions solvation phenomenon; conversely, the first one is almost unaffected. The peak positions in the g _OO(r) of pure water and supercooled-liquid state are in close agreement with the values reported in the literature [¹⁸]. On the other hand, all oxygen atoms accept two hydrogen bonds from neighboring water molecules, and nearly every hydrogen atom interacts with another oxygen atom. The coordination number of O - O for water is about 6.39 [³], four molecules form the tetrahedral arrangement via the hydrogen bonds and the two additional non-hydrogen-bonded water molecules participate in the first shell of solvate on of Cl^- o rLi⁺[³⁹].

In the solution, the first peaks of g _OH(r) and g _HH(r), centered on r = 0.97 Å (Fig. 2b) and r = 1.5 Å (Fig. 2c), respectively, are logically associated with the intramolecular OH-bond and H - H correlation. These findings are in close agreement with the experimental values and also with those calculated by simulation methods for the intern structure of pure water [¹³,¹⁸,²¹]. There is no discrepancy in the intramolecular structure of water concerning pure water, suggesting that the internal structure of the water molecule remains unchanged in any case. Therefore, neither changes of state nor the presence of ions influence the well-known structure of the water molecule.

3.2.2. Solvation structure of the ion pair

In this part, our discussions will be focused on the structure of the lithium and chlorine hydration-shells, and their pair distribution functions are respectively observed in Fig. 3a), b) and Fig. 3c), d).

Figure 3 Lithium-Oxygen a) lithium-hydrogen, b) chloride-oxygen, c) chloride-hydrogen, and d) pair distribution functions and running integration numbers, from HRMC simulation of LiCl6H₂O at the supercooled-liquid state. 

For a better understanding of various aspects of the structural behavior of aqueous solutions on the atomic or molecular level, it is very important to investigate the hydration shells. In electrolyte solutions, the polar nature of the water molecule generates particular distributions around ions to form the hydration shells. The negative partial charge of the oxygen atom is attracted by the cation Li⁺, and the positive partial charges of hydrogen atoms are repulsed by the same cation; the positively charged cations orient water molecules so that they situate non-bonded oxygen atoms near the cation [³]. In aqueous solutions, water molecules - by their large dielectric constant- can greatly reduce the electrostatic interaction between cations and anions; thus the association and dissociation of the LiCl ion pairs become easy and frequent [⁵].

a. Hydration Shell of Li⁺: The pair distribution functions of g _OLi(r) and g _HLi(r) together with the corresponding running integration numbers, obtained from our RMC [³] and HRMC simulations are shown in Fig. 3a) and b). According to the numbers and the average positions of hydrogen and oxygen first near-neighbors of lithium r = Å, the number of water molecules forming the first hydration shell of Li⁺ is about three to four water molecules [³]. The coordination number obtained here is in good agreement with those reported in the literature [¹³].

On the other hand, the two hydration shells are distinct. The first minima which follow the two first peaks of the pair distribution functions of Li-O and Li-H are well-defined; this means that the lithium first and second solvation shells are clearly separated. It is also interestingly to note that the first hydration shell is more structured compared to the second one; this is explained by the fact that the first peaks for the two pair distribution functions are quite pronounced. We can also note that our findings confirm one of results cited in literature [³⁸], there is almost no exchange between the first and second hydration shell, meanwhile the minimum that separates first and second peaks of the Li-O correlations tends to zero. It is also noteworthy that the contact ion pair formation leaves the hydration shell of Li steady and unchanged.

b. Hydration Shell of Cl^-: Now, we turn to the discussion of the structure of the hydration shells around the chlorine ion. As in the previous subsection, the pair distribution functions of g _OCl(r) and g _HCl(r), together with the corresponding running integration numbers, obtained from our RMC [³] and HRMC simulations, are represented in Fig. 3c) and d). As shown in Fig. 3c), the first peak of g _OCl(r) appears at 3.2 Å, and the behavior of running coordination number plot for g _OCl(r) indicates the existence of a defined first hydration shell at this position. These findings are confirmed by the behavior of the running coordination number plot for g _HCl(r) (see Fig. 3d), which shows a similar behavior at the same average position. It should be also noted that remarkable pre-peaks appear at the neighborhood of r = 2.2 Å for both pair distribution functions of O-Cl and H-Cl. The number of water molecules due to these pre-peaks is about 1.5 to 2. We suggest that it represents the two non-bonded hydrogen water molecules discussed in the section of the water structure.

Unlike the Li’s hydration shells, the contact ion pair formation disturbs the hydration shell of Cl by the penetration of the Li-ion with the majority of the water molecules forming its hydration shells, in the hydration shell of the chlorine.

This disturbance is even more important when the concentration is higher. We can conclude that this disturbance greatly reduced the number of water molecules in the first hydration layer, which appeared as pre-peaks in the Cl-H and Li-H distribution functions. This poor solvation is compensated by the hydration shell observed around the average position of 3.2 Å. Finally, the chlorine ion owns, as the lithium one, two hydration shells but very close and unsteady.

3.2.3. Pair distribution functions: Ion-Ion

Concerning the ion-ion correlations, no structure is observed for the same ions correlations (see Fig. 4a) and c)) because of the strong repulsions. The only significant function is that of LiCl shown in Fig. 4b). The very sharp, intense peak at r = 2.92 Å is due to the contact ion pair formation Li-Cl discussed above.

The superposition of the Cl-H and Cl-O distribution functions on that of Cl-li shows the penetration of the lithium hydration layers and further confirms the results discussed in the hydration layer part.

Figure 4 Pair distribution functions Ion-Ion: g _ClCl(r) (a), g _ClLi(r) (b) and g _LiLi(r) (c) at the supercooled-liquid state. 

4 Conclusion

We have investigated the solvation structure of the LiCl in aqueous solution using Reverse Monte Carlo simulation at 162 K. The findings emphasize the utility of potential on the Reverse Monte Carlo simulation and suggest that the supercooled-liquid solution behaves like a distinguished state compared to the others thermodynamic states liquid and glass studied in previous works. The behavioral differences for the glassy state-observed through the structural anomalies are confirmed here. At this state, the two ions of the LiCl pair own two hydration shells, but those of lithium are distinct, better defined, more steady, and represent at the same time the main reason for the unsteady of those of chlorine. The Li’s hydration shells overlap with those of chlorine and disrupt its solvation in the contact ion pair formation.