Using LAMMPS to shed light on Haven’s ratio: Calculation of Haven’s ratio in alkali silicate glasses using molecular dynamics

Haven and Verkerk studied the diffusion of ions in ionic conductive glasses with and without an external electric field to better understand the mechanisms behind ionic conductivity. In their work, they introduced the concept now known as Haven’s ratio (HR), which is defined as the ratio of the tracer diffusion coefficient (Dself) of ions to the diffusion coefficient from steady-state ionic conductivity (Dσ), calculated by the Nernst–Einstein equation. Dσ can be challenging to obtain experimentally because the number of charge carriers has to be implied, a subject still under discussion in the literature. Molecular dynamics (MD) allows for direct measurement of the mean squared displacement (r 2) of diffusing cations, which can be used to calculate D, avoiding the definition of a charge carrier. Using MD, the authors have calculated the r 2 of three alkali ions (Li, Na, and K) at different temperatures and concentrations in silicate glass, with and without the influence of an electric field. Results found for HR generally fell close to 0.6 at lower concentrations (x = 0.1) and close to 0.3 at higher concentrations (x = 0.2 and 0.3), comparable to the literature, implying that the electric field introduces new mechanisms for the diffusion of ions and that MD can be a powerful tool to study ionic diffusion in glasses under external electric fields.


Introduction
In a liquid, particles are free to flow. In contrast, a solid forms because the energy of the particles decreases, allowing them to take on a relatively ordered, three-dimensional, and more rigid structure. Even so, particles can flow within a solid, given the right circumstances. For example, in glasses, monovalent cations have been found to diffuse throughout their matrix relatively quickly. Some glass compositions have already been studied for different applications because of their remarkably high alkali ion diffusion, especially under an external electric field (Daiko et al., 2022). The interest in this particular situation goes beyond scientific curiosity. These materials may be applied in solid-state batteries and sensors, and the phenomenon is closely related to the mixed-alkali effect (MAE) (Calahoo et al., 2020).
In the case of a single charge carrier, the diffusion under the influence of an external electric field at a temperature T can be given by the ionic conductivity (σ) using the Nernst-Einstein equation (Varshneya and Mauro, 2019): where Z is the valence of the charge carrier (in the case of a monovalent cation, Z 1), e is the elementary charge, n is the density of effective charge carriers, D is the diffusion coefficient, k B is Boltzmann's constant, and T is the absolute temperature. For a given temperature and chemical composition, σ is entirely dependent on n and D. However, determining which parameters control the phenomenon can be complex, as obtaining the terms empirically and independently from each other is challenging with the currently available technology. Consequently, no universally accepted model explains the mechanisms behind the ionic conductivity in all glass materials (Dyre et al., 2009). Generally, the models found in the literature can be divided into strong electrolytes or weak electrolytes, depending on whether the ions are entirely mobile in the glass matrix or partially immobilized (Bragatto, 2020). This distinction implies that n might be equal to the total density of the species responsible for the ionic conductivity or just a fraction of the total density. Therefore, an essential step in developing such a universal model is to gain a better understanding of the nature of a charge carrier.
One way to look at this problem is by analyzing Haven's ratio (H R ) (Varshneya and Mauro, 2019). H R is defined as the ratio of the self-diffusion coefficient of ions (D Self ) to the diffusion coefficient calculated from conductivity (D σ ) using Eq. 1: ( 2 ) According to Haven and Verkerk (1965), H R can be interpreted as follows: for an ionic crystal in which the diffusion is interstitial, the values for D Self and D σ are the same, and H R = 1. On the other hand, if the mobile ion follows a different random walk, in which some of the jumps are not allowed during D Self but are allowed during D σ , the diffusion coefficient will be different, with D σ > D Self . A more indepth analysis of H R along with its meaning and impact on the mechanisms of ionic conduction in glasses can be found in the research by Isard (1999), Kahnt (1996), andMurch (1982).
In their work, Haven and Verkek found the values of H R were between 0.4 and 0.6 for different sodium silicate glasses. The authors also observed that H R decreased sharply with increasing composition up to 10 mol% and kept dropping more slowly at higher concentrations. Since their work, many authors have investigated H R , obtaining similar values for other ionicconductive glasses.
D Self can be measured in a laboratory using radioactive tracers. However, calculating D σ from Eq. 1 might be problematic as the true nature of n must be known. In general, n is taken as the total ion concentration per volume, implying that the charge carriers act as a strong electrolyte. Molecular dynamics (MD) is a novel and powerful tool that can shed some light on this problem. Previous studies (Welch et al., 2019;Atila et al., 2020;Zhao et al., 2020) were able to measure the diffusion of different species of glassy materials using MD without making assumptions about the nature of the charge carriers. In other words, D can be understood from the mean squared displacement (r 2 ) of mobile ions by where n d is the dimensionality within which the process occurs and Δt is the observational time. In turn, r 2 is given by where N is the number of averaged particles, x i (0) is an initial reference position, and x i (t) is the position with respect to time. Although experiments on the diffusion of species in glasses are common, the study of the same diffusion with an external electric field is not. In this work, glasses of the composition xA 2 O·(1-x)SiO 2 , with x = 0.1, 0.2, and 0.3, and A = Li, Na, and K were prepared using MD, and the r 2 values of the alkali ions with and without an external field were calculated. These values were used to obtain D Self , D σ , and H R . We hypothesize that studying these parameters with MD and comparing them to laboratory experimental results and interpretations in the literature will validate the use of this technique and expand the current application of MD in the glass science field to study the material under the effects of an external electric field.

Methods
All simulations were run using the Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS), an opensource molecular dynamics software application (Du and Cormack, 2022). Alkali silicate glasses of the form xA 2 O·(1-x)

FIGURE 1
Example of the resulting glass made using the methodology described in this work, with 5100 atoms of lithium (yellow), oxygen (red), and silicon (black) at the right proportions to give a composition of 0.30 Li 2 O 0.70 SiO 2 .
Frontiers in Materials frontiersin.org 02 SiO 2 were simulated for x = 0.1, 0.2, and 0.3 and A = Na, Li, and K. A visual representation can be seen in Figure 1.
Pairwise potentials developed by Pedone et al. (2006) were used to describe the interactions between the species according to where D ij is the bond dissociation energy, a ij is a function of the slope of the potential energy well, and r 0 is the equilibrium bond distance. The parameters derived by Pedone et al. from binary oxide crystals pertinent to the systems in this work are given in Table 1. These potentials were chosen due to their widely reported accuracy for silicates (Pedone et al., 2008;Pedone, 2009;Jabraoui et al., 2016). Simulation volumes were chosen to match experimentally obtained density values for each composition. The density values used to find the volume of the simulation cell for each glass composition from this work can be found in Table 2 18 . A time step of 1 femtosecond (fs) was used for each composition. All simulations were split into two steps. The first step was the creation of the simulated glass, and the second step included testing each composition at various temperatures with and without an applied electric field. A benefit of splitting the simulations into steps was that an almost identical glass structure was present for each simulation. Each composition was simulated independently in phase one, with 5100 atoms correctly proportioned. Atoms were placed randomly inside the simulation cell, with densities, shown in Table 2, determining the volume.
In all simulations, a microcanonical ensemble that constrained the number of particles, volume, and internal energy (NVE) was used to facilitate energy distribution in the simulation cell. After initialization with an NVE ensemble for 100 ps, the simulations were introduced to a temperature of 3000 K in an NVT ensemble (constant number of particles, volume, and temperature) for 100 ps. Next, with energy equilibrium being reached, the pressure was allowed to equilibrate through the NPT ensemble (constant number of particles, pressure, and temperature) for 100 ps. The glasses were then cooled to 500 K at a rate of 1 K/ps. Here, each glass composition was duplicated among the various temperature and electric field tests.
Phase two started with another NVE initialization process for 100 ps. Then, the glass was raised to T at 1 K/ps in an NPT ensemble  Pedone et al. (2006).

FIGURE 2
Schematic representation of preparation steps for measurement at different temperatures. Different temperatures were used to produce a reliable glass structure starting from a random distribution of the atoms. Beyond the vertical line, an electric field of 0.05 V/Å was applied to the x-axis.

FIGURE 3
r 2 as a function of time for the glass 0.3 Na 2 O 0.7 SiO 2 at 1000 K. In this simulated experiment, a field of 0.05 V/A was applied in the x-axis. The dotted lines represent the best linear fit (R 2 ≥ 0.997), as found by the OriginPro 2019 software. The first 100 ps were discarded for this calculation to guarantee a steady-state diffusion regime.
Frontiers in Materials frontiersin.org 03 where T = 800 K, 1000 K, 1500 K, 2000 K, 2500 K, and 3000 K. Finally, an external electric field of 0.05 V/Å was applied across the x dimension of the simulation cell. This field strength was chosen as lower values had no detectable difference between exposed and unexposed samples. A schematic representation of the experiment is shown in Figure 2. In comparison, higher values resulted in unrealistic values for MSD and diffusion of the less mobile atoms, possibly resulting from dielectric breakdown. Similar electric field strengths were observed in other simulation works (Welch et al., 2019).
Another valid and more common concern in this kind of work is the time length of the diffusional process (Du and Cormack, 2022). Due to molecular dynamic limitations, the observational time is very short compared to laboratory experiments. Therefore, in our calculations of diffusion coefficient, we discarded the first 100 ps of every run to avoid any non-diffusive behavior by the ions.

Results
For each glass composition and each temperature, six different experiments were completed, measuring r 2 as a function of time. Three were completed without any electric field applied, and three experiments applied an electric field of 0.05 V/Å on the x-axis. An example can be seen in Figure 3, with results for the glass 0.3 Na 2 O 0.7 SiO 2 at 1000 K.
To obtain reliable values for the diffusion coefficient D, the diffusion must be at a steady-state regime (Varshneya and Mauro, 2019), meaning that the dependency of r 2 with time should be linear. To determine if the regime could be considered at a steady-state pace, all data were fitted using OriginPro 2019 software, and all results had R 2 ≥ 0.997 after disregarding data points obtained before the 100 ps. Values of r 2 Δt were obtained from these fittings, and values for D Self and D σ were calculated using Eq. 3. Also, it was possible to obtain values of H R using Eq. 2. An example of these results for measurements determined at 1000 K can be found in Table 3.
Values for diffusion coefficients found in this work were comparable to laboratory diffusion coefficients for sodium in soda-lime silicate glasses (Mehrer et al., 2008), previous simulation results found by Welch et al. (2019), and thermal diffusion values found using MD in the literature (Du and Chen, 2012;Du and Cormack, 2022). This should be taken with caution because an artifact of the molecular dynamics experiment is that the glass's temperature and the fictive temperature might be underestimated (Du and Cormack, 2022).
Values for activation energy (E A ) can be obtained from the dependency of the diffusion coefficient with temperature, assuming an Arrhenius behavior (Varshneya and Mauro, 2019): where D can be either D Self or D σ , and D 0 is a constant. Results for this work can be found in Figure 4, expressed as log D as a function of 1000/T. The activation energies were obtained by fitting the data in Figure 4 and can be found in Table 4. Activation energy values for the self-diffusion were comparable to values found in the literature for lithium silicate glasses using MD (Du and Chen, 2012).

Discussion
The data obtained to calculate Haven's ratio for all glass compositions at the temperature interval range described earlier are shown in Figure 4. The results can be divided into two categories: composition dependence and temperature dependence.

Compositional dependence
The most significant differences between the three alkali-oxide glasses were their ionic radii, electronegativity, and binding energies. As shown in Table 4, K showed a large discrepancy compared to Li and Na. To illustrate this, Figure 5 shows the data for diffusion coefficients at 1000 K. This result might seem contradictory because Li is known for being the best alkali conductor of the three alkalis considered in this research (Otto and Milberg, 1968). However, it presented a diffusion coefficient smaller than K and similar to Na through the compositional range.  3 Values for D Self and D σ obtained using linear fits from r 2 results and Eq. 3 for measurements determined at 1000 K. The errors were obtained by the standard deviation of three different results for D self and three different results for D σ from different MD experiments.

FIGURE 4
Logarithm of the diffusion coefficient as a function of the inverse of temperature for the different glasses studied in this work. The linear fitting (R 2 ≥ 0.95) indicated an Arrhenius behavior and was used to obtain the activation energies for the different measurements. The error was calculated from the linear regression (OriginPro 2019).
TABLE 4 Activation energies obtained using linear fits from the average values for D Self and D σ shown in Figure 4 (R 2 ≥ 0.95) and Eq. 6. Errors presented were obtained from the mathematical fit (OriginPro 2019).

Glass modifier
Glass modifier concentration Another exciting aspect that can be seen in Figure 5 was the similarity between the results found for Li and Na-containing glasses. This result disagreed with experimental results, in which the difference in conductivity was about one order of magnitude (Otto and Milberg, 1968). Ionic conductivity can be interpreted by Eq. 1 and simplified as the ionic conductivity of a given glass composition is proportional to the product of the density of effective charge carriers (n) and the diffusion coefficient (D). In the experiments, the diffusion coefficients were calculated, but the number of effective charge carriers was not. This was due to the nature of the simulation experiment, where the short observed time frame was comparable to voltage frequencies of 10 9 Hz. In comparison, the laboratory ionic conductivity results were obtained with voltage frequencies ranging from 10 -2 to 10 7 Hz (Irvine et al., 1990).
Following this logic, the similarity between the molecular dynamic diffusion coefficients and the difference between their laboratory ionic conductivity might be justified as a difference between the number of effective charge carriers. Furthermore, this approach suggests a "weak-electrolyte" behavior by the glass, in which the dissociation equilibrium plays a more significant role in the conductivity. More on the two theories can be found elsewhere (Martin and Angell, 1986).

FIGURE 5
Self-diffusion and electric field induced diffusion coefficients for the different glasses included in this research at 1000 K. Results are comparable to values found in the literature (Mehrer et al., 2008;Welch et al., 2019).

FIGURE 6
Haven's ratio (H R ) for the different glasses used in this research at 1000 K. Values of H R are higher for low concentrations (x = 0.1) and decrease with increasing alkali concentration (x = 0.2 and 0.3). Similar behavior was found for other glass systems (Isard, 1999) and crystals (Murch, 1982) in laboratory experiments.

Frontiers in Materials
frontiersin.org 06 Independent of the reasons for differences in diffusion coefficients, H R was calculated using the diffusion coefficients obtained, as described in the Methods section. For illustration, the results obtained at 1000 K are presented in Figure 6. As can be seen, the values for low (x = 0.1) alkali concentrations were close to 0.6, and for higher concentrations (x = 0.2 and 0.3), the values were found to be between 0.1 and 0.3. This fast decrease of H R with composition is known in the literature for different oxide glass systems (Kelly et al., 1980;Thomas and Peterson, 1984;Bychkov et al., 2001) and has two possible explanations: either the diffusional mechanisms change from low concentrations to high concentrations, or the mechanisms are the same, but the ion-ion or defect-ion interactions are different (Bychkov et al., 2001).
To the best of the authors' knowledge, there is no definitive experimental confirmation for either the weak versus strong electrolyte approach or the change of diffusional mechanisms versus the ion-ion or ion-defect interaction approach. At the same time, MD studies could help shed light on these issues and provide new information that would be impossible to obtain otherwise.

Temperature dependence
The variation of the diffusion coefficients with temperature can be seen in Figure 4. The difference between the values was more significant at lower temperatures and decreased with the increase in temperature, especially at higher alkali concentrations. This behavior was observed for all the glasses studied in this work, where the two diffusion coefficients were equal at higher temperatures for the Li and Na glasses and close to equal for the K glasses. The values were reasonable, considering that the activation energies for D self are higher than D σ , but should be equal when the system is in the liquid state and the ionic charge carrier is fully dissociated (Garrido et al., 2018). The difference between D self and D σ was more easily expressed by H R , as described in Eq. 2.

FIGURE 7
Values of Haven's ratio found in this work using MD for the glasses x A 2 O (1-x) SIO 2 , with x = 0.1, 0.2, or 0.3, and A = Li, Na, or K.
Frontiers in Materials frontiersin.org 07 Results for H R as a function of temperature for all the glasses studied in this work are found in Figure 7. Overall, the absolute value of H R increased with the temperature, but at a lower rate for the low concentration (x = 0.1) and a faster pace for the higher concentrations (x = 0.2 and 0.3).
The same small positive temperature dependence was observed by Isard and was attributed to a distribution of activation energies among the vacancies that the alkali used to diffuse (Isard, 1999). A similar dependence on temperature was observed in crystalline systems (Murch, 1982). This was not the only interpretation of the dependency of H R with temperature found in the literature, as some authors presented H R to be an indication of preferred pathways (Bychkov et al., 2001), different mechanisms for the diffusion with or without an external electric field (Haven and Verkerk, 1965), or even mechanisms of correlated motion (Kahnt, 1996). Similar to its dependency on composition, the authors believe that MD experiments can be a powerful tool to help researchers better understand H R and its meaning.
From the good agreement of both compositional and temperature dependency of Haven's ratio on the studied glass system, the present work strongly indicated that molecular dynamics is a powerful tool to better understand this phenomenon. Many of the interpretations in the literature (Murch, 1982;Kahnt, 1996;Isard, 1999) associated the changes in the transport mechanisms with structural variations in the vicinity of the mobile ions and, consequently, a distribution of activation energies. Both structure and activation energies can be obtained using MD and will be presented and discussed by the authors in future works.

Conclusion
Haven's ratio and its origin are not yet fully understood, and obtaining good and reliable results can be challenging as it requires assumptions about the nature of ionic conductivity and expensive measurements for self-diffusion. However, as a relatively new technique, molecular dynamics can be a fast, cheap, and easy addition to research, given the proper considerations. In this work, the authors have examined simple alkali silicate glasses to verify the validity of MD techniques in this context. Due to its nature, MD simulation experiments have limitations, especially with sample size and observational times, both relevant to the ionic conductivity phenomenon. Therefore, to observe a response from the system, the temperatures and voltage field applied had to be comparably higher than in laboratory experiments. The shift to higher temperatures is well-known in the literature. Still, the higher electric field used in this work was found by trial and error, and it was enough to create a significant diffusion of the alkali studied but not affect the diffusion of the glass-former units.
Values of r 2 were obtained using the LAMMPS software, and from them, values for D self and D σ were calculated. Their dependency on temperature, alkali concentration, and alkali nature was also studied. Results obtained in this work agree with different points found in the literature, such as the compositional dependency of H R , where the value sharply decreases at concentrations above x = 0.1 of alkali oxide, and temperature dependence, in which H R approaches unity at higher temperatures. These results are encouraging and show the possibility of using MD techniques to further investigate the diffusional properties of oxide glasses, especially when under an external electric field.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.