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

Abstract

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 (H_R), which is defined as the ratio of the tracer diffusion coefficient (D_self) 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²) 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² 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 H_R 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.

1 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 is the valence of the charge carrier (in the case of a monovalent cation, , is the elementary charge, n is the density of effective charge carriers, is the diffusion coefficient, is Boltzmann’s constant, and is the absolute temperature.

For a given temperature and chemical composition, is entirely dependent on and . 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 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 () (Varshneya and Mauro, 2019). is defined as the ratio of the self-diffusion coefficient of ions () to the diffusion coefficient calculated from conductivity () using Eq. 1:

According to Haven and Verkerk (1965), can be interpreted as follows: for an ionic crystal in which the diffusion is interstitial, the values for and are the same, and = 1. On the other hand, if the mobile ion follows a different random walk, in which some of the jumps are not allowed during but are allowed during , the diffusion coefficient will be different, with > . A more in-depth analysis of 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), and Murch (1982).

In their work, Haven and Verkek found the values of were between 0.4 and 0.6 for different sodium silicate glasses. The authors also observed that decreased sharply with increasing composition up to 10 mol% and kept dropping more slowly at higher concentrations. Since their work, many authors have investigated , obtaining similar values for other ionic-conductive glasses.

can be measured in a laboratory using radioactive tracers. However, calculating from Eq. 1 might be problematic as the true nature of must be known. In general, 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, can be understood from the mean squared displacement () of mobile ions bywhere is the dimensionality within which the process occurs and is the observational time. In turn, is given bywhere is the number of averaged particles, is an initial reference position, and 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₂O⋅(1-x)SiO₂, with x = 0.1, 0.2, and 0.3, and A = Li, Na, and K were prepared using MD, and the values of the alkali ions with and without an external field were calculated. These values were used to obtain , and . 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.

2 Methods

All simulations were run using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), an open-source molecular dynamics software application (Du and Cormack, 2022). Alkali silicate glasses of the form xA₂O∙(1-x)SiO₂ 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.

FIGURE 1

Pairwise potentials developed by Pedone et al. (2006) were used to describe the interactions between the species according towhere D_ij is the bond dissociation energy, a_ij is a function of the slope of the potential energy well, and r₀ 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¹⁸.

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 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).

FIGURE 2

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.

3 Results

For each glass composition and each temperature, six different experiments were completed, measuring 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₂O 0.7 SiO₂ at 1000 K.

FIGURE 3

To obtain reliable values for the diffusion coefficient , the diffusion must be at a steady-state regime (Varshneya and Mauro, 2019), meaning that the dependency of 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² ≥ 0.997 after disregarding data points obtained before the 100 ps. Values of were obtained from these fittings, and values for and were calculated using Eq. 3. Also, it was possible to obtain values of 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 () can be obtained from the dependency of the diffusion coefficient with temperature, assuming an Arrhenius behavior (Varshneya and Mauro, 2019):where can be either or , and is a constant. Results for this work can be found in Figure 4, expressed as as a function of . The activation energies were obtained by fitting the data in Figure 4 and can be found in Table 4.

FIGURE 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).

4 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.

4.1 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.

FIGURE 5

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 (). 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⁹ Hz. In comparison, the laboratory ionic conductivity results were obtained with voltage frequencies ranging from 10^–2 to 10⁷ 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).

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).

FIGURE 6

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.

4.2 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.

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).

FIGURE 7

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.

5 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² 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.

References

• 1

  AtilaA.OuaskitS.HasnaouiA. (2020). Ionic self-diffusion and the glass transition anomaly in aluminosilicates. Phys. Chem. Chem. Phys.22, 17205–17212. 10.1039/d0cp02910f

  • CrossRef
  • Google Scholar

• 2

  BansalN. P.DoremusR. H. (2013). Handbook of glass properties. Elsevier Science.

  • Google Scholar

• 3

  BragattoC. B. (2020). “Electrical properties of glass,” in The world scientific reference of amorphous materials (World Scientific), Vol. 15, 199–222.

  • Google Scholar

• 4

  BychkovE.PriceD. L.LappA. (2001). Universal trend of the haven ratio in glasses: Origin and structural evidences from neutron diffraction and small-angle neutron scattering. J. Non-Crystalline Solids293–295, 211–219. 10.1016/s0022-3093(01)00673-1

  • CrossRef
  • Google Scholar

• 5

  CalahooC.XiaY.ZhouR. (2020). Influence of glass network ionicity on the mixed‐alkali effect. Int. J. Appl. Glass Sci.11, 396–414. 10.1111/ijag.15546

  • CrossRef
  • Google Scholar

• 6

  DaikoY.SakudaA.HonmaT.HayashiA. (2022). Ionic conduction of glasses and their potential applications. J. Ceram. Soc. Jpn.130, 552–557. 10.2109/jcersj2.22035

  • CrossRef
  • Google Scholar

• 7

  DuJ.ChenC.-H. (2012). Structure and lithium ion diffusion in lithium silicate glasses and at their interfaces with lithium lanthanum titanate crystals. J. Non.-Crys. Sol.358, 3531–3538. 10.1016/j.jnoncrysol.2012.05.015

  • CrossRef
  • Google Scholar

• 8

  DuJ.CormackA. N. (2022). Atomistic simulations of glasses. John Wiley and Sons. 10.1002/9781118939079.fmatter

  • CrossRef
  • Google Scholar

• 9

  DyreJ. C.MaassP.RolingB.SidebottomD. L. (2009). Fundamental questions relating to ion conduction in disordered solids. Rep. Prog. Phys.72, 046501. 10.1088/0034-4885/72/4/046501

  • CrossRef
  • Google Scholar

• 10

  GarridoL.AranazI.GallardoA.GarciaC.GarciaN.BenitoE.et al (2018). Ionic conductivity, diffusion coefficients, and degree of dissociation in lithium electrolytes, ionic liquids, and hydrogel polyelectrolytes. J. Phys. Chem. B122, 8301–8308. 10.1021/acs.jpcb.8b06424

  • CrossRef
  • Google Scholar

• 11

  HavenY.VerkerkB. (1965). Phys. Chem. Glasses6, 38.

  • Google Scholar

• 12

  IrvineJ. T. S.SinclairD. C.WestA. R. (1990). Electroceramics: Characterization by impedance spectroscopy. Adv. Mater.2, 132–138. 10.1002/adma.19900020304

  • CrossRef
  • Google Scholar

• 13

  IsardJ. O. (1999). The Haven ratio in glasses. J. Non-Crystalline Solids246, 16–26. 10.1016/s0022-3093(99)00036-8

  • CrossRef
  • Google Scholar

• 14

  JabraouiH.VaillsY.HasnaouiA.BadawiM.OuaskitS. (2016). Effect of sodium oxide modifier on structural and elastic properties of silicate glass. J. Phys. Chem. B120, 13193–13205. 10.1021/acs.jpcb.6b09664

  • CrossRef
  • Google Scholar

• 15

  KahntH. (1996). Ionic transport in glasses. J. Non-Crystalline Solids203, 225–231. 10.1016/0022-3093(96)00354-7

  • CrossRef
  • Google Scholar

• 16

  KellyJ. E.CordaroJ. F.TomozawaM. (1980). Correlation effects on alkali ion diffusion in binary alkali oxide glasses. J. Non-Crystalline Solids41, 47–55. 10.1016/0022-3093(80)90190-8

  • CrossRef
  • Google Scholar

• 17

  MartinS. W.AngellC. A. (1986). Dc and ac conductivity in wide composition range Li2O-P2O5 glasses. J. Non-Crystalline Solids83, 185–207. 10.1016/0022-3093(86)90067-0

  • CrossRef
  • Google Scholar

• 18

  MehrerH.ImreÁ.Tanguep-NijokepE. (2008). Diffusion and ionic conduction in oxide glasses. J. Phys. Conf. Ser.106, 012001. 10.1088/1742-6596/106/1/012001

  • CrossRef
  • Google Scholar

• 19

  MurchG. E. (1982). The haven ratio in fast ionic conductors. Solid State Ionics7, 177–198. 10.1016/0167-2738(82)90050-9

  • CrossRef
  • Google Scholar

• 20

  OttoK.MilbergM. E. (1968). Ionic conduction in alkali and thallium silicate glasses. J. Am. Ceram. Soc.51, 326–329. 10.1111/j.1151-2916.1968.tb15947.x

  • CrossRef
  • Google Scholar

• 21

  PedoneA.MalavasiG.MenzianiM. C.CormackA. N.SegreU. (2006). A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses. J. Phys. Chem. B110, 11780–11795. 10.1021/jp0611018

  • CrossRef
  • Google Scholar

• 22

  PedoneA.MalavasiG.MenzianiM. C.SegreU.CormackA. N. (2008). Molecular dynamics studies of Stress−Strain behavior of silica glass under a tensile load. Chem. Mat.20, 4356–4366. 10.1021/cm800413v

  • CrossRef
  • Google Scholar

• 23

  PedoneA. (2009). Properties calculations of silica-based glasses by atomistic simulations techniques: A review. J. Phys. Chem. C113, 20773–20784. 10.1021/jp9071263

  • CrossRef
  • Google Scholar

• 24

  ThomasM. P.PetersonN. L. (1984). Electrical conductivity and tracer diffusion in sodium germanate glasses. Solid State Ionics14, 297–307. 10.1016/0167-2738(84)90114-0

  • CrossRef
  • Google Scholar

• 25

  VarshneyaA. K.MauroJ. C. (2019). Fundamentals of inorganic glasses. Elsevier.

  • Google Scholar

• 26

  WelchR. S.WilkinsonC. J.MauroJ. C.BragattoC. B. (2019). Charge carrier mobility of alkali silicate glasses calculated by molecular dynamics. Front. Mat.6, 121. 10.3389/fmats.2019.00121

  • CrossRef
  • Google Scholar

• 27

  ZhaoY.DuJ.QiaoX.CaoX.ZhangC.XuG.et al (2020). Ionic self-diffusion of Na2O–Al2O3–SiO2 glasses from molecular dynamics simulations. J. Non-Crystalline Solids527, 119734. 10.1016/j.jnoncrysol.2019.119734

  • CrossRef
  • Google Scholar