# 1 Introduction

## 1 Introduction

In recent years several molecular dynamics computer simulations have been done in order to investigate the structure and dynamics of sodium silicate melts and glasses (Smith, Greaves, and Gillan 1995, Cormack and Cao 1997). By using the potential proposed by Vessal, Amini, Fincham and Catlow (1989) these authors have found that the structure of systems like, e.g., sodium disilicate (SDS) is characterized by a microsegregation in which the sodium atoms form clusters of a few atoms between bridged SiO<sub>4</sub> units. In order to see whether this somewhat surprising result is reproduced also by a different model from the one of Vessal et al. (1989) we have performed simulations of SDS using a different potential (discussed in detail below). In addition to the investigation of the structure we also study the dynamical properties of SDS in order to see whether the finite size effects that have been observed in pure silica (Horbach, Kob, Binder, and Angell 1996, Horbach, Kob, and Binder 1999a, and Horbach, Kob, and Binder 1999b) are present in SDS as well.

## 2 Model

The model potential we use to describe the interactions between the ions in SDS is the one proposed by Kramer, de Man, and van Santen (1991) which is a generalization of the so called BKS potential (van Beest, Kramer, and van Santen 1990) for pure silica. It has the following functional form:

$$\varphi (r)=\frac{q_\alpha q_\beta e^2}{r}+A_{\alpha \beta }\mathrm{e}^{B_{\alpha \beta }r}\frac{C_{\alpha \beta }}{r^6}\alpha ,\beta [\mathrm{Si},\mathrm{Na},\mathrm{O}].$$
(1)   
Here $`r`$ is the distance between an ion of type $`\alpha `$ and an ion of type $`\beta `$. The values of the parameters $`A_{\alpha \beta },B_{\alpha \beta }`$ and $`C_{\alpha \beta }`$ can be found in the original publication. The potential (1) has been optimized by Kramer et al. for zeolites, i.e. for systems that have Al ions in addition to Si, Na and O. In that paper the authors used for silicon and oxygen the partial charges $`q_{\mathrm{Si}}=2.4`$ and $`q_\mathrm{O}=1.2`$, respectively, whereas sodium was assigned its real ion charge $`q_{\mathrm{Na}}=1.0`$. With this choice charge neutrality is not fulfilled in systems like SDS. To overcome this problem we introduced for the sodium ions a position dependent charge $`q(r)`$ instead of $`q_{\mathrm{Na}}`$,

$$q(r)=\{\begin{array}{cc}0.6\left(1+\mathrm{ln}\left[C\left(r_\mathrm{c}r\right)^2+1\right]\right)\hfill & r<r_\mathrm{c}\hfill \\ 0.6\hfill & rr_\mathrm{c}\hfill \end{array}$$
(2)   
which means that for $`rr_\mathrm{c}`$ charge neutrality is valid ($`q(r)=0.6`$ for $`rr_\mathrm{c}`$). Note that $`q(r)`$ is continuous at $`r_\mathrm{c}`$. We have fixed the parameters $`r_\mathrm{c}`$ and $`C`$ such that the experimental mass density of SDS and the static structure factor from a neutron scattering experiment (see below) are reproduced well. From this fitting we have obtained the values $`r_\mathrm{c}=4.9`$Å and $`C=0.0926`$Å<sup>-2</sup>. With this choice the charge $`q(r)`$ crosses smoothly over from $`q(r)=1.0`$ at $`1.7`$ Å to $`q(r)=0.6`$ for $`rr_\mathrm{c}`$.

The simulations have been done at constant volume with the density of the system fixed to $`2.37\mathrm{g}/\mathrm{cm}^3`$. The equations of motion were integrated with the velocity form of the Verlet algorithm and the Coulombic contributions to the potential and the forces were calculated via Ewald summation. The time step of the integration was $`1.6`$ fs. In this paper we investigate the properties of SDS in the liquid state at $`T=2100`$ K and in the glass state at $`T=300`$ K. The equilibration time at $`T=2100`$ K was two million time steps thus corresponding to a real time of $`3.5`$ ns. At this temperature we simulated systems with $`N=1008`$ and $`N=8064`$ particles. In order to improve the statistics two independent runs were done for the large system and eight independent runs for the small system. The glass state was produced by cooling the system from equilibrated configurations at $`T=1900`$ K with a cooling rate of $`1.1610^{12}\mathrm{K}/\mathrm{s}`$. The pressure is $`4.5`$ GPa at $`T=2100`$ K and $`0.96`$ GPa at $`T=300`$ K.

## 3 Results

In order to demonstrate that our model is able to reproduce the structure of real SDS very well we compare the static structure factor $`S^{\mathrm{neu}}(q)`$ with the one from a neutron scattering experiment by Misawa, Price, and Suzuki (1980). To calculate $`S^{\mathrm{neu}}(q)`$ one has to weight the partial structure factors from the simulation with the experimental coherent neutron scattering lengths $`b_\alpha `$ ($`\alpha [\mathrm{Si},\mathrm{Na},\mathrm{O}]`$):

$$S^{\mathrm{neu}}(q)=\frac{1}{_\alpha N_\alpha b_\alpha ^2}\underset{kl}{}b_kb_l\mathrm{e}^{i𝐪[𝐫_k𝐫_l]}.$$
(3)   
The values for $`b_\alpha `$ are $`0.414910^{12}`$ cm, $`0.36310^{12}`$ cm and $`0.580310^{12}`$ cm for silicon, sodium and oxygen, respectively. They are taken from Susman, Volin, Montague, and Price (1991) for silicon and oxygen and from Bacon (1972) for sodium. Fig. 1 shows $`S^{\mathrm{neu}}(q)`$ from the simulation and the experiment at $`T=300`$ K. We see that the overall agreement between simulation and experiment is good. For $`q>2.3`$ Å<sup>-1</sup>, which corresponds to length scales of next nearest Si–O and Na–O neighbors, the largest discrepancy is at the peak located at $`q=2.8`$ Å<sup>-1</sup> where the simulation underestimates the experiment by approximately $`15`$% in amplitude. Very well reproduced is the peak at $`q=1.7`$ Å<sup>-1</sup>, which is called the first sharp diffraction peak, and which is a prominent feature in pure silica as well. In silica this peak arises from the tetrahedral network structure since the length scale which corresponds to it, i.e. $`2\pi /1.7\mathrm{\AA }^1=3.7`$ Å, is approximately the spatial extent of two connected SiO<sub>4</sub> tetrahedra. From the figure we recognize that this structure is partly present in SDS also. The peak at $`q=0.95`$ Å<sup>-1</sup> is not present in the experimental data which might be due to the fact that in this $`q`$ region the experimental resolution is not sufficient. By looking at the coordination number distributions, discussed below, we see that the peak at $`q=0.95`$ Å<sup>-1</sup> is related to a super structure which is formed by the sodium and silicon atoms. In agreement with this interpretation the length scale corresponding to this peak, i.e. $`2\pi /0.95\mathrm{\AA }^1=6.6`$ Å, is two times the mean distance of nearest Na–Na or Na–Si neighbors.

The coordination number distribution $`P_{\alpha \beta }(z)`$ for different pairs $`\alpha \beta `$ gives the probability that an ion of type $`\alpha `$ has exactly $`z`$ nearest neighbors of type $`\beta `$. By definition two neighboring atoms have a distance from each other which is less than the location of the first minimum $`r_{\mathrm{min}}`$ of the corresponding partial pair correlation function $`g_{\alpha \beta }(r)`$. From the functions $`g_{\alpha \beta }(r)`$ we find for $`r_{\mathrm{min}}`$ the values $`3.6`$ Å, $`5.0`$ Å, $`2.35`$ Å, $`5.0`$ Å, $`3.1`$ Å, and $`3.15`$ Å for the Si–Si, Si–Na, Si–O, Na–Na, Na–O, and O–O correlations. We note that $`P_{\mathrm{Si}\mathrm{O}}(z)`$ is larger than $`0.99`$ for $`z=4`$ at $`T=2100`$ K and $`T=300`$ K which means that nearly every silicon atom is four fold coordinated with oxygen atoms forming a SiO<sub>4</sub> tetrahedron. Some of the distribution functions are shown in Fig. 2. We recognize from Fig. 2a that about $`65`$% of the oxygen atoms are bridging oxygens between two tetrahedra ($`P_{\mathrm{O}\mathrm{Si}}(z=2)0.65`$), and that about $`28`$% of the oxygen atoms form dangling bonds ($`z=1`$) with corresponding silicon atoms. In the neighborhood of these dangling bonds sodium atoms are located. This means that the sodium atoms partly destroy the (disordered) tetrahedral network which is the structure for pure silica. A significant number of oxygen atoms have no silicon atoms, but only sodium atoms as direct neighbors. In Fig. 2b we show $`P_{\mathrm{Na}\mathrm{Na}}(z)`$ at $`T=2100`$ K and $`T=300`$ K, and we see that essentially the two distributions coincide with a a mean value between $`z=8`$ and $`z=9`$. Basically the same distribution is found for $`P_{\mathrm{Na}\mathrm{Si}}(z)`$. Therefore, every sodium atom is surrounded by $`8`$$`9`$ other sodium atoms and $`8`$$`9`$ silicon atoms. Since the mean distance between Na–Na and Na–Si neighbors is approximately the same, i.e. about $`3.3`$ Å, we can conclude that sodium and silicon atoms form a spherical super structure in which every sodium atom is surrounded by a first shell of oxygen atoms and a second shell of silicon and sodium atoms. In Fig. 2b we have also included $`P_{\mathrm{Na}\mathrm{O}}(z)`$ and we recognize that every sodium atom has on average about $`4`$$`5`$ nearest oxygen neighbors. There is again no essential difference between $`P_{\mathrm{Na}\mathrm{O}}`$ for $`T=2100`$ K and $`P_{\mathrm{Na}\mathrm{O}}`$ for $`T=300`$ K, although the pressure is about a factor $`4.5`$ higher at $`T=2100`$ K. This means that the structural features which we observe at $`T=2100`$ K are not formed due to the relatively high pressure. Nevertheless, we emphasize that the small difference between $`P(z)`$ in the liquid and in the glass state is partly due to the high cooling rate of about $`10^{12}\mathrm{K}/\mathrm{s}`$ we used to produce the structures at $`T=300`$ K. A careful analysis of the cooling rate effects for SDS will be presented elsewhere (Horbach, Kob, and Binder 1999c).

Having described the structure of SDS we turn now our attention to a dynamical quantity, namely the self part of the dynamic structure factor $`S_\mathrm{s}(q,\nu )`$, which depends on frequency $`\nu `$ and the magnitude of the wave–vector $`q`$. It is defined by

$`S_\mathrm{s}(q,\nu )`$ $`=`$ $`{\displaystyle \frac{N_\alpha }{N}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t\mathrm{e}^{2\pi \nu t}\mathrm{e}^{i𝐪[𝐫_\alpha (t)𝐫_\alpha (0)]}`$ (4)   
$`=`$ $`{\displaystyle \frac{N_\alpha }{N}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t\mathrm{e}^{2\pi \nu t}F_\mathrm{s}(q,t)\alpha [\mathrm{Si},\mathrm{Na},\mathrm{O}]`$   
where $`𝐫_\alpha (t)`$ is the position vector of a tagged particle of type $`\alpha `$ at time $`t`$ and $`F_\mathrm{s}(q,t)`$ is the incoherent intermediate scattering function. The details of the Fourier transformation in (4) are given elsewhere (Horbach et al. 1999a). $`S_\mathrm{s}(q,\nu )`$ for silicon, sodium and oxygen is shown in Fig. 3a at the temperature $`T=2100`$ K and the wave–vector $`q=1.7`$ Å<sup>-1</sup> for the two system sizes $`N=1008`$ and $`N=8064`$. For sodium the curves for the small and the large system coincide over the whole frequency range. This is not the case for silicon and oxygen for which $`S_\mathrm{s}(q,\nu )`$ has no system size dependence for frequencies $`\nu 1.1`$ THz whereas for smaller frequencies there is a missing of intensity for the small system. The vibrational modes causing a shoulder around $`\nu =0.9`$ THz in $`S_\mathrm{s}(q,\nu )`$ are usually called boson peak excitations. In the small system a part of these excitations is missing due to the loss of intensity for $`\nu 1.1`$ THz. An explanation of this behavior for the case of silica, which shows qualitatively the same finite size effects, can be found in Horbach et al. (1999a and 1999b). Due to the sum rule $`𝑑\nu S_\mathrm{s}(q,\nu )=1`$ the missing of the boson peak excitations for frequencies 0.4 THz $`\nu `$ 1.1 THz in the small system has to be “reshuffled” to smaller frequencies leading to a broadening and an increase of the quasielastic line around $`\nu =0`$. Since the quasielastic line is outside the frequency resolution of our Fourier transformation the consequences in the change of the quasielastic line can be observed better in the Fourier transform of $`S_\mathrm{s}(q,\nu )`$, i.e. the incoherent intermediate scattering function $`F_\mathrm{s}(q,t)`$, which is shown in Fig. 3b for the system sizes $`N=8064`$ and $`N=1008`$ at $`q=1.7`$ Å<sup>-1</sup>. We recognize from this figure that $`F_\mathrm{s}(q,t)`$ shows a two step relaxation behavior similar to the case of silica and fragile glassformers (Ngai, Riande, and Ingram 1998). As expected from our results for $`S_\mathrm{s}(q,\nu )`$ the scattering functions $`F_\mathrm{s}(q,t)`$ have no system size dependence for sodium. In $`F_\mathrm{s}(q,t)`$ for silicon and oxygen the height of the plateau increases and the $`\alpha `$ relaxation process shifts to longer times with decreasing system size. Furthermore, the scattering functions for the small system show a pronounced oscillation for $`t>0.2`$ ps which is due to the fact that the boson peak excitations present in the small system cause a peak in $`S_\mathrm{s}(q,\nu )`$ whereas in the large system only a shoulder is present. Finally we mention that the finite size effects in the dynamics of SDS are found in the whole $`q`$ range and, moreover, do not affect the static properties of SDS.

Acknowledgments:

This work was supported by SFB 262/D1 and by Deutsch–Israelisches Projekt No. 352–101. We also thank the RUS in Stuttgart for a generous grant of computer time on the Cray T3E.

## References

* Bacon, G. E., 1972, Acta Cryst. A, 28, 357.
* Cormack, A. N., and Cao, Y., in Modelling of Minerals and Silicated Materials, Eds.: B. Silvi and P. Arco (Kluwer, Dordrecht 1997).
* Horbach, J., Kob, W., Binder, K., and Angell C. A., 1996, Phys. Rev. E, 54, R5889.
* Horbach, J., Kob, W., and Binder, K., 1999a (submitted to Phys. Rev. B)
* Horbach, J., Kob, W., and Binder, K., 1999b, to appear in the Proceedings on Neutrons and Numerical Methods, Grenoble, 1998, preprint in cond-mat/9901162
* Horbach, J., Kob, W., and Binder, K., 1999c (to be published)
* Kramer, G. J., de Man, A. J. M., and van Santen, R. A., 1991, J. Am. Chem. Soc., 113, 6435.
* Misawa, M., Price, D. L., and Suzuki, K., 1980, J. Non–Cryst. Solids, 37, 85.
* Ngai, K. L., Riande, E., and Ingram, M.D., (editors), 1998, J. Non–Cryst. Solids, 235–237 (Proceedings of the Third International Discussion Meeting on Relaxations in Complex Systems, Vigo, 1997).
* Smith, W., Greaves, G. N., and Gillan, M. J., 1995, J. Chem. Phys., 103, 3091.
* Susman, S., Volin, K. J., Montague, D. G., and Price, D. L., 1991, Phys. Rev. B, 43, 11076.
* van Beest B. W., Kramer G. J., and van Santen R. A., 1990, Phys. Rev. Lett., 64 1955.
* Vessal, B., Amini, A., Fincham, D., and Catlow, C. R. A., 1989, Philos. Mag. B, 60, 753.