Solvent Exchange around Aqueous Zn(II) from Ab Initio Molecular Dynamics Simulations

: Hydrated zinc(II) cations, due to their importance in biological systems, are the subject of ongoing research concerning their hydration shell structure and dynamics. Here, ab initio molecular dynamics (AIMD) simulations are used to study solvent exchange events around aqueous Zn 2+ , for which observation in detail is possible owing to the considerable length of the generated trajectory. While the hexacoordinated Zn(H 2 O) 62+ is the dominant form of Zn(II) in an aqueous solution, there is a non-negligible contribution of the pentacoordinated Zn(H 2 O) 52+ complex which presence is linked to the dissociative solvent exchange events around Zn 2+ . The pentacoordinated Zn(II) has a much tighter hydration sphere and is characterized by a trigonal bipyramidal structure, in contrast to the usual octahedral symmetry of the hexacoordinated complex. In total, two full exchange events are registered in the analyzed trajectory. AIMD simulations on an adequate length scale thus provide a direct way of studying such solvent exchange events around ions in molecular detail.


Introduction
The coordination of zinc(II) cation is of paramount importance in biological systems, as zinc is the active center in hundreds of metalloenzymes, most often catalyzing hydrolysis reactions [1]. The importance of zinc goes beyond the role of cofactor in proteins, as the zinc cation is also actively utilized for intracellular signaling and subject to elaborate homeostasis mechanisms in the cell [2].
Due to the particular importance of water as the solvent in biological context, the structure and dynamics of hydrated ions is a topic of special significance in physical chemistry of solutions [3][4][5]. The influence of ions on water is often expressed in the terms of structural effects exerted by ions on the solvent, leading to the distinction between structure making and breaking ions [5].
AIMD based on DFT is a particularly suited technique for studying aqueous ionic solutions [51,52]. Beginning with the pioneering study of Be 2+ (aq) [53], many metal ions in an aqueous solution have been studied using AIMD simulations [54][55][56][57][58]. The current availability of fast and approximate DFT methods, such as density functional tight binding (DFTB) allows for the generation of extensive data sets on solvated ions in molecular solvents, including water [59]. The particular advantage of AIMD is the ability to study intermolecular polarization effects that ultimately influence many dynamic properties of the studied systems [54,55,60].
Computational methods are also indispensable for elucidating solvent exchange mechanisms around ions [61,62]. Experimental methods, rather than providing a mechanism of the exchange reaction, can only help in estimating the rate of this process [24]. The details of H 2 O exchange around Zn(II) have been studied using MD simulations and static QM calculations [28][29][30]37]. Although AIMD simulations, combining the quantum description of electrons with the classical propagation of nuclei according to Newton's laws of motion, seem ideal for this purpose, the rather long timescale of water exchange in the first hydration shell of Zn(II) has thus far hampered direct observation of this process using AIMD.
In this work, in order to bridge this gap and obtain reliable data on the structure of Zn 2+ (aq) and the dynamics of the solvent exchange process in its first hydration shell, we use extensive DFT-based AIMD simulations with a computational setup that is proven to be especially suitable for investigating liquid water. We hypothesize that the inability of hitherto equilibrium AIMD studies to register any solvent exchange around Zn(II) was a limitation of the short timescale of the simulations, rather than the flows in the method itself. We also aim to register the interplay of well-defined solvation complexes of Zn 2+ (aq) that participate in the solvent exchange process.

Computational Methods
We studied aqueous Zn(II) in a system composed of a single Zn 2+ cation and 100 H 2 O molecules, which was contained in a cubic simulation supercell with applied periodic boundary conditions. The volume of the system was chosen to represent the experimental density of H 2 O at 298 K (997 kg m -3 [63]) coupled with the absolute standard partial molar ionic volume of Zn 2+ (aq) (−34.9 cm 3 mol -1 [64,65]), resulting in cell volume 14.33 3 Å 3 . The starting configuration was prepared using GROMACS (v. 2018.5) tools [66] by solvating a centrally placed Zn 2+ cation with a solvent configuration obtained from a well-equilibrated AMOEBA water simulation [67].
The system was initially subject to energy minimization and subsequently equilibrated using classical MD simulations in the canonical (NVT) ensemble at T = 298.15 K performed using Tinker-HP package (v. 1.2) [68]. The revised AMOEBA force field for liquid water [69] in combination with the AMOEBA parameterization of Zn(II) [38] were applied. The ultimate configuration from this NVT trajectory was then used to initialize the actual AIMD simulation.
AIMD simulations were performed with the cp2k computational suite (v. 6.0) [70][71][72], employing the DFT-based representation of the electronic structure as implemented in the Quickstep electronic structure module in cp2k [73]. The revPBE exchange-correlation functional [74,75] coded in the libXC library [76] was used, as it provides excellent reproduction of numerous static and dynamic properties of liquid water and aqueous ionic solutions [77][78][79][80][81]. A mixed Gaussian atomic orbitals with plane waves (GPW) representation of the electronic structure is applied in cp2k [82], and we used a molecularly optimized short-ranged double-zeta (DZVP-MOLOPT-SR-GTH) basis set for atomic orbitals [83], coupled with the auxiliary plane wave expansion of the electron density up to a 500 Ry cutoff. Though heavily contracted and featuring few diffuse primitive atomic orbitals, the chosen basis set performs favorably when compared to a more expensive triple-zeta polarized basis set for the description of liquid water [84]. Only valence electrons were treated explicitly, while core electrons were represented by norm-conserving GTH pseu-dopotentials [85] parameterized for the PBE functional. Additionally, the smoothing of the electron density and its derivative on the spatial integration grid was applied (via keywords XC_SMOOTH_RHO NN50 and XC_DERIV NN50_SMOOTH in cp2k), as it was previously found to significantly improve the stability of the local energetics of liquid water [86]. Dispersion effects were included via two-body DFT-D3 empirical dispersion correction with zero damping term [87] with the cutoff set to 15 Å. The studied system has a non-zero total charge (+2) and the charge compensation is provided by an implicit neutralizing jellium background [88], as implemented in cp2k.
The system was first equilibrated for ∼35 ps in an AIMD simulation in the canonical (NVT) ensemble at T = 298.15 K using the canonical velocity rescaling (CSVR) thermostat [89] with the time constant set to 2000 cm −1 (≈16.67 fs). The simulation time step was set to 0.5 fs. After this equilibration period, the thermostat was removed and a production AIMD simulation in the microcanonical (NVE) ensemble was continued for 100 ps, with data collection every 2 fs. The centers of maximally localized Wannier functions (MLWFs) [90] were also computed and stored every 2 fs. Analyses were performed using in-house code and VMD (v. 1.9.4) [91], which also served as a visualization tool, along with gnuplot (v. 5.2) [92].

Results
The static structure of the hydration shells of aqueous Zn(II) is summarized by radial distribution functions (RDFs) for the Zn· · · O and Zn· · · H pairs, as depicted in Figure 1. The studied system is large enough to fully contain two hydration shells of Zn 2+ , as clearly seen in the Zn· · · O RDF that features a prominent shallow minimum with g(r) ≈ 0 separating the first hydration shell from the second one. At first glance, the running integration number of g ZnO (r) is very close to 6 along this minimum, implying a predominantly hexacoordinated Zn(II) ion. However, its exact value of 5.8 suggests that during as much as 20% of the trajectory Zn 2+ might in fact be pentacoordinated. Detailed parameters of the analyzed RDFs are gathered in Table 1.   Figure 1: location of the first and the second maxima and minima and the respective running integration numbers.

RDF
The first hydration shell of Zn 2+ extends up to 3.03 Å and contains a non-integer number of water oxygen atoms. In order to determine the instantaneous coordination number (CN) of Zn(II), one could simply count O atoms found closer than the r min,1 value, which would be tantamount to accepting the Heaviside step function as a selector function for H 2 O molecules. However, this results in a discontinuous CN trajectory and adopting a smooth selector function instead can help visualize transitions between different coordination states of Zn 2+ (aq), as previously proposed [50]. In this work, the coordinated solvent molecules are selected using a logistic cutoff function: controlled by the cutoff radius r 0 and the sharpness parameter s. r ZnO,i (t) is the instantaneous Zn· · · O distance for the water molecule i. Adopting r 0 = 3.1 Å and s = 0.08 Å ensures a smooth transition between a fully coordinated H 2 O molecule in the first hydration shell, where CN i (t) → 1, and a non-coordinated molecule in the second hydration shell, where CN i (t) → 0, cf. Figure 1.
The time evolution of thus obtained smooth coordination number of Zn 2+ is shown in Figure 2. As clearly seen, even though a continuous CN(t) definition is accepted, its value provides a clear-cut distinction between the hexa-and pentacoordinated states of aqueous Zn(II). The zinc cation spends most of the time in a Zn(H 2 O) 6 2+ complex (about 81.2% of the trajectory). The pentacoordinated complex may be unanimously found in ∼16.3% of the trajectory, while we ascribe the rest to transitional structures or short-time existence periods of one of the two dominating structures that are too short for inclusion in the analysis. Armed with the knowledge about the existence of two distinct coordination states of Zn 2+ (aq), we revisit Figure 1 in order to compare the Zn· · · O RDFs in the regions of existence of the hexa-and pentacoordinated complexes. As seen in Figure 3, there is only a slight difference between the RDF based on the entire trajectory and the one based only on the fragments where Zn(H 2 O) 6 2+ is found, which is unsurprising given the predominance of the hexacoordinated Zn(II). More importantly, the Zn· · · O RDF restricted to the trajectory periods where Zn(H 2 O) 5 2+ exists reveals a major compression of not only the first but also the second hydration shell of the zinc cation. The first maximum shifts to 2.07 Å, while the second one to 4.07 Å. Thus, the two distinct complexes are characterized with different hydration patterns extending up to the entire second hydration shell.  Unlike static geometry optimizations, AIMD simulations can provide ensembleaveraged structures. In order to find the preferred geometry of hexa-and pentacoordinated Zn 2+ (aq), we isolated the first shell complexes from the relevant trajectory fragments, rotated them to a common reference frame and averaged the geometry. The results can be seen in Figure 4. Zn(H 2 O) 6 2+ forms a regular octahedral complex with the average Zn· · · O distance equal 2.14 Å. On the contrary, Zn(H 2 O) 5 2+ is characterized by a trigonal bipyramidal arrangement of water ligands. The Zn· · · O distance for the three equatorial molecules is 2.04 Å, while the axial O atoms are located on average 2.14 Å from the metal cation. The geometries of the two complexes in Cartesian coordinates are available in the Supplementary Materials, see Figure S1. While the positions of O atoms tend to be fixed in the Zn(II) solvation cage, the hydrogen atoms of the first shell water molecules are relatively free to bend from the coplanar arrangement with the Zn-O vector. We define the tilt angle, β, as the angle formed between the Zn-O vector and the dipole moment vector of the coordinated H 2 O molecule. The two-dimensional distribution of water molecule configurations in the (r ZnO , β) plane is shown in Figure 5. The distribution of angles is evidently bimodal, with a global maximum at 0°, indicating a coplanar arrangement of the H 2 O molecular plane and the Zn-O vector, and a local maximum at ca. 30°, reflecting the presence of a substantial population of water molecules with tilted H atoms. We note here that the latter tend to be located slightly farther away from the central atom. The presence of two complexes in aqueous solution of Zn(II) has far-reaching implications for the water molecules in its hydration shell. Ions are known to polarize water and the extent of this polarization is related to the changes in infrared spectrum of solution with respect to bulk water [54,55]. In Figure 6 we compare the probability distribution functions of the length of molecular dipoles of H 2 O molecules in the entire system vs the ones in the first hydration shell of Zn 2+ (aq). Unlike the dipole moment in bulk water, which has a symmetric Gaussian distribution [54,55,60], the curve for the studied system is right-skewed due to the substantial polarization of the neighboring H 2 O molecules by the zinc cation. The distribution can be fitted to a split Pearson VII peak shape with a maximum located at 2.78 D. On the contrary, the distributions of dipole moments of water molecules found in the two complexes are symmetric and can similarly be fit to a Pearson VII analytical shape. It is found that H 2 O molecules in Zn(H 2 O) 5 2+ are polarized more significantly (on average 3.71 D) than the ones in Zn(H 2 O) 6 2+ (mean 3.43 D). We finally would like to estimate the time frame of a single exchange event based on the trajectories of the participating water molecules. As seen in Figure 2, we register two major exchange events connected to a change in CN of Zn(II), at around 0-12 ps and again at 78-86 ps. Additionally, a short-lived occurrence of Zn(H 2 O) 5 2+ is found at ca. 64 ps. In order to shed light on the molecular details of these exchange events, we plotted the time evolution of the Zn· · · O distance of involved H 2 O molecules, as presented in Figure 7.
The first event begins as soon as the production trajectory starts. The water molecule W 47 dissociates from the initially present Zn(H 2 O) 6 2+ at 0.72-0.75 ps. Hereafter, the W n notation indicates a H 2 O molecule with an arbitrary index n in the trajectory. The solvent exchange process is obviously dissociative, as no new solvent molecule coordinates to the cation for an extended period of time. In this event, the released H 2 O moves to the second hydration shell, as its separation from Zn 2+ increases to ca. 4 Å (cf. Figure 1). This results in a transient Zn(H 2 O) 5 2+ complex that persists up to ∼8.8 ps, at which moment an interesting concerted mechanism is observed. Namely, W 47 is kicked off to the weakly organized third hydration shell of the cation, at 6 Å and beyond, while another second shell H 2 O, specifically W 12 , coordinates to zinc(II), restoring the Zn(H 2 O) 6 2+ complex. This, however, proves short-lived, as W 12 dissociates already at ∼10.1 ps and returns to the second hydration shell. Since a similar phenomenon is observed in further events, vide infra, such attempted exchange events that ultimately fail to stabilize the newly formed complex are most probably a common occurrence in Zn 2+ (aq). Somewhat surprisingly, the entire exchange event is finalized by a third water molecule (W 22 ) that is captured by Zn 2+ not from the second, but rather directly from the third hydration shell, initially >6 Å away from the cation. The hexacoordinated complex is restored at ∼11.9 ps and remains stable for an extended period of time. Interestingly, the coordination of W 22 is seemingly again assisted by W 47 , since the time of arrival of the latter in the third hydration shell coincides with the start of movement of the former towards the inner hydration shells. The second exchange event is only an attempted escape of W 9 from Zn(H 2 O) 6 2+ at ca. 64 ps. However, no other H 2 O molecule tries to enter the first hydration shell of Zn(II) this time and ultimately, after transiently wandering off to the verge of the second shell, W 9 is recaptured by the complex and the stable octahedral arrangement of ligands persists.
The third exchange event (and the second successful one) resembles the first one to a great extent. It begins with the dissociation of W 77 at ca. 78.3 ps. This H 2 O molecule is released initially to the second hydration shell, but quickly travels further to ultimately stay at ca. 5 Å separation from zinc(II) for an extended period, which is the location of the shallow minimum between the second and the poorly defined third hydration shell of the cation (cf. Figure 1). Interestingly, this perturbation of the hydration shells leads to two concerted movements of water molecules towards the central complex, initially that of W 40 , quickly followed by W 49 . Similarly to the first exchange event, the intervening water molecule (here, W 40 ) ultimately does not manage to become properly coordinated to the cation. Instead, it remains at the very edge of the first hydration shell, only sligthly perturbing the newly formed Zn(H 2 O) 5 2+ complex. After its release back into the second hydration shell at ca. 81 ps, the pentacoordinated complex remains stable for almost 4 ps. Then, an interesting cascade of concerted events can be noticed. W 77 returns to the second hydration shell, apparently swapping its place with W 40 , but at the same time W 49 manages to get captured by zinc(II), reinstating the hexacoordinated complex. Again, concerted movement of H 2 O molecules seems to be crucial for the mutually assisted exchanges between solvation shells leading to coordination number change of Zn 2+ (aq).
The particular advantage of our analysis is the rigorous separation of the trajectory fragments with hexa-vs. pentacoordinated Zn(II), thanks to the smooth instantaneous coordination number, see Figure 2. This allowed us to prove the tightening of the hydration shells of Zn 2+ Figure 3) are in good agreement with earlier AIMD simulations, a similar separation of trajectory fragments was proposed (2.11 Å and 2.05 Å, respectively [50]). Static geometry optimizations also lead to a general conclusion that the water molecules are located closer to Zn(II) in Zn(H 2 O) 5 2+ than in Zn(H 2 O) 6 2+ [26,27,29,31,32]. The hexacoordinated complex always maintains the octahedral arrangement of water oxygens (although the inclusion of hydrogens lowers the point group to T h [9]), which is also found in our time-averaged structure, cf. Figure 4. More controversial is the preferred geometry of Zn(H 2 O) 5 2+ . Based on geometry optimizations, both trigonal bipyramid (point group D 3h ) [26,31,32] and square pyramid (point group C 2v ) [27,29,94] have been confirmed as minima of the potential energy surface. It appears though that the presence of second shell H 2 O molecules stabilizes the former structure [29,40,94]. The presence of external hydration shells is thus crucial to the observation of the time-averaged trigonal bipyramidal structure that we report here, see Figure 4. The discrepancy between the Zn· · · O distance to equatorial and axial molecules (2.04 Å and 2.14 Å, respectively) is similar to that found previously in an isolated cluster (2.05 Å and 2.13 Å, respectively [26]).
Notably for AIMD simulations, our considerable system size with 100 H 2 O molecules allows us to record the RDF for the entire second hydration shell, as well as for the weakly formed third shell, see Figure 1. Due to the structure-making nature of Zn 2+ (aq), the second shell is also accessible experimentally. Previous EXAFS results indicate a second shell of 12.2 O atoms at an average distance of 4.40 Å [16] or 11.2 O atoms at 4.10 Å [14]. The distance to second shell from XRD measurements equals 4.02-4.26 Å [20][21][22]. Results from AIMD and QM/MM MD simulations generally support these experimental findings, with second hydration shell composed of 13-15 molecules at an average distance 4.3-4.5 Å [42,43,48,50]. In comparison, our global average of the position of the second peak in g ZnO (r) (4.11 Å, see Figure 1 and Table 1) is generally in the lower limit of the reported distances, while the number of O atoms (15.2) indicates a fairly crowded second hydration shell of Zn(II). Most importantly, the lowering of CN of the cation results in a major compression of this shell, with the RDF peak shifted to 4.07 Å (see Figure 3) and at least two water molecules released to the bulk solution, as the number of O atoms in the second shell decreases to 13.
Earlier studies on monovalent ions found that they are strongly polarized by water, but in turn do not polarize the solvent significantly (in the sense of only slight alteration of H 2 O dipole moment) [54,55,60]. In strike contrast to these findings, the water molecules in the first hydration shell of Zn 2+ (aq) are heavily polarized by the electric field of the cation, see Figure 6. For all H 2 O molecules, we find an average dipole moment of 2.78 D, which is slightly lower than the value usually found in AIMD simulations of liquid water (∼3 D) [54,60]. However, mean dipole moment of liquid H 2 O for revPBE/revPBE0 functionals is indeed slightly lower and in agreement with the present findings [96,97]. While the increase of dipole moment of H 2 O near the cation was previously found in MD simulations with the polarizable Amoeba force field, as well as AIMD simulations [38,48], the clear advantage of our analysis is the ability to separate the two forms of hydrated Zn(II). In the hexacoordinated complex, the dipole moment of the first shell water molecules increases by 0.65 D with respect to the bulk. Dissociation of this complex leads to a further increase of the average dipole moment and it is almost 1 D greater than for non-coordinated H 2 O. Thus, the water molecules in Zn(H 2 O) 5 2+ are evidently more heavily polarized than in Zn(H 2 O) 6 2+ . The 0.6 D increase of the mean dipole moment of first hydration shell H 2 O has been noted previously in AIMD simulations [48]. Simulations with the Amoeba force field indicate that the water molecules closest to the cation have their dipole moments raised by more than 1 D [38].
Computational methods are excellent tools for studying solvent exchange mechanism around solvated ions [62]. The experimental estimate of water residence time on Zn(II) from quasi-elastic neutron scattering is a rather broad range of 0.1-5 ns [24]. Assuming the actual residence time is near the lower limit, it is in principle accessible to longer MD simulations. Previously, only force field description of the potential energy surface could provide the possibilities to reach the required time scale. However, the exact formulation of Zn(II)-H 2 O interactions is crucial for enabling solvent exchange events and they were frequently not observed even in prolonged MD simulations [16,[34][35][36]. On the other hand, a water residence time of 2.2 ns was found in an MD simulation with the Amoeba force field [38] and a rate constant for water exchange equal 4.1 × 10 8 s -1 was also obtained [37]. Most AIMD or QM/MM MD simulations to date suffered from length limitations and no exchange events were observed [42][43][44][45][46][47][48]. In stark contrast to these earlier results, we managed to register two full solvent exchange events in the first hydration shell of Zn(II). While it is impossible to quantitatively estimate the water residence time from such a small sample, we note that the starts of successful exchange events are separated by ca. 78 ps, so that our results are not in disagreement with the experimental estimates provided the true residence time is close to the lower limit proposed previously [24]. We note here that the solvent exchange may be studied using accelerated sampling, e.g., metadynamics [37,50,93]; however, while such an approach provides the details of the potential energy surface of the hydrated cation, only equilibrium simulations can capture the actual timescale of the observed events.
Unfortunately, the question of how the two complexes differ on the free energy scale cannot be answered directly from our equilibrium simulation. However, we can roughly estimate this free energy difference from the fractions of the trajectory attributed to the penta-and hexacoordinated structures as ∆G = −RT ln(0.163/0.812) ≈ 4 kJ mol -1 . On the one hand, this is high enough that the pentacoordinated structure might prove too elusive to be registered experimentally, e.g., with EXAFS, noting also that the average coordination number of Zn(II) (5.8) is within the typical uncertainty of structure determination methods from the value of 6.0, indicative of perfect hexacoordination. On the other hand, this is low enough (only 1.6RT) to allow relatively free interconversion between the two structures. We note here that the free energy difference between [UO 2 (H 2 O) 5 ] 2+ (aq) and [UO 2 (H 2 O) 4 ] 2+ (aq) from ab initio metadynamics is ∼3 kJ mol -1 , with similarly inconclusive experimental results as in the case of Zn(II) [98].
Our results indicate that the solvent exchange around Zn(II) follows a dissociative mechanism, cf. Figure 7. The release of a single water molecule from the first hydration shell leads to a transient Zn(H 2 O) 5 2+ complex, the lifetime of which may be tentatively estimated as 8-10 ps, based on the limited sample of events we registered. A similar timescale of dissociative exchange events was noted before by Fatmi et al. [43]. A similar mechanism of solvent exchange around hexacoordinated Zn(II) was also inferred from static QM calculations [29,30,62]. Interestingly, the dissociation of a water molecule does not lead to increased CN in the second shell. On the contrary, it undergoes compression and major structural changes, as evident from Figure 3. Consequently, the registered exchange events require concerted movement of assisting solvent molecules, cf. Figure 7.
This cooperative movement of H 2 O molecules enables the capture of a newly coordinated water molecule from as far as the third solvation shell of Zn(II).

Conclusions
In this paper, the structure of the hydrated zinc(II) cation has been studied using extensive AIMD simulations. In contrast to numerous earlier MD investigations at different levels of theory, which frequently found only hexacoordinated zinc with fixed first hydration shell waters, we clearly identify the presence of both hexa-and pentacoordinated complexes of Zn 2+ . While the octahedral Zn(H 2 O) 6 2+ remains the dominant form of Zn(II) in an aqueous solution, Zn(H 2 O) 5 2+ is present in about 16% of the trajectory. The average structure of the latter is a trigonal bipyramid and its presence results in major structural changes reaching up to the second hydration shell and beyond. Radial distribution functions indicate major compression of the second hydration shell of the pentacoordinated complex and molecular dipole distributions show profound polarization of first shell H 2 O molecules that is also heavily dependent on CN.
Owing to the considerable length of the generated trajectory, we were able to register two full solvent exchange events around aqueous Zn 2+ . They undoubtedly follow a dissociative mechanism and are driven by concerted solvent movement that may feature molecules up to the third hydration shell of Zn(II). While, due to the limited sampling, we cannot quantitatively estimate the water residence time on the zinc cation, the separation of successful events by ca. 80 ps seems to be consistent with the lower limit of experimental residence times (about 0.1 ns). The single exchange event, accompanied by the presence of Zn(H 2 O) 5 2+ , seems to take about 8-10 ps. We conclude that equilibrium AIMD simulations on an adequate length scale provide a direct way of studying solvent exchange events around aqueous ions in molecular detail. The smooth selector function for coordinated solvent is a powerful tool that allows the identification of solvent exchange events and their preferred (associative or dissociative) mechanism. The analysis framework presented here may be tremendously useful for other solvated ions and more complicated systems containing Zn(II), such as concentrated zinc halide solutions.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
Calculations were performed at the Academic Computer Center in Gdańsk (TASK).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: AIMD ab initio molecular dynamics CN coordination number DFT density functional theory DFTB density functional tight binding EXAFS extended X-ray absorption fine structure GTH Goedecker-Teter-Hutter (pseudopotentials) MD molecular dynamics MLWF maximally localized Wannier function ND neutron diffraction QM quantum mechanics QM/MM quantum mechanics/molecular mechanics revPBE revised Perdew-Burke-Ernzerhof (functional) XRD X-ray diffraction