Local microstructure of silica glass

Local microstructure of silica glass

Journal of Non-Crystalline Solids 355 (2009) 1215–1220 Contents lists available at ScienceDirect Journal of Non-Crystalline Solids journal homepage:...

620KB Sizes 0 Downloads 9 Views

Journal of Non-Crystalline Solids 355 (2009) 1215–1220

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol

Local microstructure of silica glass L.T. Vinh a, P.K. Hung a,*, N.V. Hong a, T.T. Tu b a b

Department of Computational Physics, Hanoi University of Technology, Viet Nam High Performing Computing Center, Hanoi University of Technology, Viet Nam

a r t i c l e

i n f o

Article history: Received 18 March 2007 Received in revised form 19 April 2009 Available online 13 June 2009 PACS: 61.43.Fs 61.20.Lc 64.70.Pf

a b s t r a c t We present a microstructure study of different amorphous states of SiO2, which are constructed by compress–decompress procedure. The local microstructure is analyzed through the coordination number, bond-angle statistic and the characteristic of void. The simulation reveals that the number of SiO4, SiO5 and SiO6 units varies in different silica states, but their topology is identical. The densest model (highest density) is obtained upon compress pressure around 20 GPa. The change in void space for different silica states is also calculated and discussed here. Ó 2009 Elsevier B.V. All rights reserved.

Keywords: Glass transition Pressure effects Microstructure Porosity Molecular dynamics Silica Medium-range order Short-range order

1. Introduction Silica is of great interest from both experimental and theoretical points of view, because it is a major component in many applications, such as semiconductor device, optical fiber and certain adsorption materials [1,2]. Recently, a numerous study focuses on so-called polyamorphism of silica. This phenomenon refers to the occurrence of distinct amorphous forms of the same substance and is observed in various systems such as GaAs, SiO2, H2O and GeO2 [3,4]. It is well known that silica glass has a tetrahedral network-structure at ambient pressure. Applying high pressure to the glass causes a formation of six-coordinated group SiO6, as is found in a high-pressure crystalline phase. In volume–pressure plots, the glass initially undergoes a volume decrease analogous to the lowdensity crystal. Within a narrow pressure range, the volume drops rapidly to value like one of the high-density crystalline solid [5–8]. Furthermore, the origin density is not fully recovered event the glass has been stored at ambient pressure for many years [9]. This experimental result implies that silica glass can exist at various stable polymorphs like those observed in crystalline silica [10].

* Corresponding author. Tel.: +84913035541. E-mail address: [email protected] (P.K. Hung). 0022-3093/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2009.05.012

Computer simulation is proven to be useful tool to study microstructure of glass and liquid [16–21]. In accordance to simulation in Ref. [11] low-density amorphous (LDA) phase SiO2 is composed primary of tetrahedral coordinated silicon and high-density amorphous (HDA) phase is 19% denser than LDA phase (in compared to 20% experimentally [9,12]). The transformation from tetrahedral at LDA to octahedral network structure at HDA is evidenced in many MD simulations [13–15,19,21]. It is found that a significant structural modification occurred at pressure window between 3 and 5 GPa is accompanied with rebounding and relaxation [15]. Furthermore, the amorphous–amorphous transition could be reversible or irreversible depending on concrete combination of pressure and temperature at which the transition takes place [13,14]. Such transition occurs gradually, but its behavior is unlike ones observed in [22]. The change in coordination number upon densification of silica glass is discussed in many works, but the systematic analysis of two terms: basic unit SiOx and void for different silica states have not been done to date. Therefore, the present paper focuses on this problem. The paper is organized as follows: in Section 2 we provide a method of preparing a silica model and of calculation of void. The result of microstructure and void in different states of silica is presented in Section 3. We end with a number of conclusions in Section 4.

1216

L.T. Vinh et al. / Journal of Non-Crystalline Solids 355 (2009) 1215–1220

2. Calculation method MD simulations of SiO2 model containing 1998 atoms (666 Si and 1332 O) in the simple cube with periodic boundary conditions are carried out. For SiO2 system a number of inter-atomic potential have been developed [11,13,19,23,24] and we adopt BKS potential, because it is still simple and well reproduces a number of physical properties in both crystal and liquid SiO2. The detail about BKS potential can be found in [24]. The long-range Coulomb interactions are calculated with the standard Ewald technique. MD time step is equal to 0.46 fs. The initial configuration is generated by randomly placing all atoms in the simulation box. This initial model is heated up to 5000 K and then treated over 50 000 steps. After that, the sample is cooled down to 3000 K and relaxed within 50 000 MD steps. Finally, the obtained liquid is cooled down again to 600 K and relaxed at constant temperature and pressure (i.e. NPT simulation) over 100 000 MD steps. This relaxation is enough long to obtain a well equilibrated model (model M1). In order to evaluate the size effect we also construct a model of 3000 atoms. No considerable difference in structural characteristics of two models has been detected. Next, we perform a compress–decompress procedure. At compress stage the model M1 is compressed to pressures Pomp upon constant temperature. Then the compressed model is relaxed at constant pressure and constant temperature for 50 000 MD steps. In second stage the model obtained at first stage is decompressed up to ambient pressure and relaxed again over 500 000 MD steps. In order to ensure that the properties of relaxed model is unchanged with further relaxation, its structural characteristics have been determined many times during last 10 000 MD steps. The positional and angular characteristics of constructed models are

(a)

calculated by averaging over 1000 configurations separated by 10 MD steps. If every atom is considered as a sphere, then there is a part of the model in which no atomic sphere lies. The radius of Si and O atoms are adopted to be 1.46 and 0.73 Å, respectively. The void is defined as a sphere that is in contact with four atoms and does not intersect with any atom. The schematic illustration of void calculation is shown in Fig. 1. The void is computed as follows: firstly, we took all set of four atoms with condition that the distance between any two atoms among them is less than 9.5 Å. Then a void sphere is inserted in contact with given atomic spheres. If the inserted void overlaps with any atom, it will be removed from system. Thereby, we obtain a set of voids, which are not intersected with any atom. Among obtained voids, some voids may be located almost inside other void or several nearest voids. These voids do not play any considerable role and therefore, next step is to remove them from system. To remove that voids, two million points are randomly generated in the simulation cell. For every ith void, we determined the number of points ni, located inside ith void. Then the number of point n2i, being in both ith void and other void is calculated. The ith void is removed if the relation n2i/ni greater than 0.95. The detail about the method of void calculation can be found elsewhere [33,34]. Void is distributed not separately, but they could group into a large cluster consisting of two or more voids. Two kinds of void aggregations are examined here: the void cluster (VC) and void tube (VT) (see Fig. 1(b) and (c)). First one is a set of voids including a central void and other voids, which overlap with central void and smaller than central void. The second one contains a number of voids bigger than oxygen atom. Furthermore, every void in VT must overlap at least with one adjacent void by a section circle

(b) Atom Atom

Void

Atom

Void

Void

Void

Void

Void Atom Atom

(d)

(c) Void

Void

Void

Central void

Void Void Void

Void

Void

Void

Void

Void

(e)

(f)

O

O

O

O Si

O

O Si

Si

Si

O O

O

O

O

O O

Fig. 1. Schematic illustration of void calculation : (a) a void is in contact with four atoms and not overlapped with any atom; (b) the voids located almost in large void (left) or in two nearest voids (right) will removed from system; (c) void cluster (VC); (d) void tube (VT); (e) one-oxygen connectivity; (f) two-oxygen connectivity.

1217

L.T. Vinh et al. / Journal of Non-Crystalline Solids 355 (2009) 1215–1220

Pressure (GPa)

5

2

70 60

4

50

3

40

Si-Si Pcomp=0 GPa Pcomp=20 GPa Pcomp=65 GPa

2

30

3

1

1

20 10

0 20

4

Si-O

0

4

5

6

7

8

Pcomp=0 GPa Pcomp=20 GPa Pcomp=65 GPa

15

3

gij(r)

Volume (cm /mol) Fig. 2. Pressure versus volume of silica glass. Branches 1, 2 correspond to compress stage when glass is compressed to 65 GPa (1) and relaxed upon this pressure (2). Branches 3, 4 related to decompressing glass to ambient pressure and its relaxation to well-equilibrium model.

10 5 0

with radius bigger than radius of oxygen atom. Such, voids of VT link to each other and create a long tube along which oxygen can move without intersection with any atom. For amorphous solids it is commonly accepted that if void radius is close to atomic radius, this void can be considered as vacancy [25–27]. Accordingly, binary amorphous system must have two kinds of vacancy, corresponding to component sizes. In the current work we examine the oxygen vacancy-like (OVV) and silicon vacancy-like void (SVV). The OVV and SVV are void with radius bigger than radius of oxygen and silicon, respectively.

O-O

4

Pcomp=0 GPa Pcomp=20 GPa Pcomp=65 GPa

3 2 1 0 0

2

4

6

8

10

-1

r(10 nm) Fig. 4. The partial radial distribution functions of silica glass.

3. Result and discussion

3.2

3.1. Atomic correlation

3

Density (g/cm )

3.1

Fig. 2 shows the pressure–volume plot as model M1 is compressed to 65 GPa and then decompressed to ambient pressure. It can be seen that the path followed by the glass is irreversible and the origin density is not fully recovered, but it approached to higher value. Similar trend is observed for another compress pressure Pcomp. Using such compress–decompress procedure we prepare six models Mx, where x = 2, 3. . .7 corresponding to Pcomp = 5, 10, 15, 20, 40 and 65 GPa. The density of model Mx versus Pcomp is displayed in Fig. 3. A maximum is seen around a point of 3.15 g/cm3. Most dense model M5 is 17% denser than model M1. Hereafter, we assume that model M5 and M1 correspond to HDA and LDA state, respectively. The structural characteristics of con-

3.0 2.9 2.8 2.7 2.6 0

10

20

30

40

50

60

70

Pcomp (GPa) Fig. 3. Density versus compress pressure for silica glass.

Table 1 Structural characteristics of silica glass. rij, gij – the position and height of first peak in PRDF; Zij – the averaged coordination number. Here 1–1 is for the Si–Si pair; 1–2 for the Si–O pair; 2–1 for the O–Si pair; 2–2 for the O–O pair. The height gij and position rij vary in the range of gij ± Dgij and rij ± Drij, where Dgij are 0.082, 0.053 and 0.046 for Si–Si, Si–O and O–O pairs, respectively; Drij is ±0.02 Å. Model

Pcomp GPa

rij, Å

gij

Zij

1–1

1–2

2–2

1–1

1–2

2–2

1–1

1–2

2–1

2–2

M1 M2 M3 M4 M5 M6 M7

– 5 10 15 20 40 65

3.12 3.12 3.10 3.12 3.12 3.10 3.10

1.60 1.60 1.60 1.60 1.60 1.60 1.60

2.58 2.58 2.56 2.54 2.54 2.56 2.56

4.00 3.89 3.60 3.39 3.38 3.57 3.52

20.02 19.13 15.57 13.24 12.89 15.59 15.95

4.10 4.01 3.53 3.29 3.28 3.47 3.47

4.33 4.41 5.00 5.52 5.58 5.04 4.96

4.07 4.10 4.28 4.43 4.45 4.25 4.24

2.03 2.05 2.14 2.22 2.23 2.13 2.12

10.06 10.23 11.31 12.00 12.11 11.28 11.11

Refs. [28,29] Ref. [23]

– –

3.12 3.15

1.62 1.59

2.65 2.59

– –

– –

– –

– –

4.00 –

2.00 –

– –

1218

L.T. Vinh et al. / Journal of Non-Crystalline Solids 355 (2009) 1215–1220 Table 2 The oxygen connectivity distribution. m is the number of common oxygen atoms that two neighbor units SiOx are bonded to. The next columns indicate the percentage of connectivity. For example, 2.499% of connectivity is two-oxygen for model 2. The statistic error is ±2.26%.

1.0

Fraction

0.8 0.6

SiO4 SiO5 SiO6

0.4 0.2

10

20

30

40

50

60

Model M1

M2

M3

M4

M5

M6

M7

1 2 3

98.431 1.546 0.023

97.494 2.499 0.007

93.616 6.285 0.099

90.711 9.034 0.255

90.013 9.747 0.240

94.813 5.093 0.094

94.871 5.056 0.073

tutes the network structure. The second angle provides the connection between two adjacent units SiOx and therefore it describes the medium range order. Fig. 6 displays O–Si–O and Si– O–Si BAD. We can see a pronounced peak appeared at 105 ± 3o for O–Si–O BAD of unit SiO4. This data is close to experimental value reported in [29] and indicates a slightly distorted tetrahedron with a silicon atom at the center and four oxygen atoms at the vertices. For units SiO5 and SiO6 there are two peaks: the main peak is located at 90 ± 3o and a small one at 140 ± 3o. It is interesting to note that the O–Si–O BADs are almost unchanged for different models. This result differs from those reported in [19], where the O–Si–O BAD gradually changes with density. The reason is that we separately calculate O–Si–O BAD for every kind among SiO4, SiO5 and SiO6 and therefore, we obtain three O–Si–O BADs instead one as provided in [19]. Such, the O–Si–O BAD of unit SiOx and averaged bond-length Si–O of them are unchanged for different models. This means that seven considered models corresponding to different states of silica glass have a network structure with identical building blocks SiO4, SiO5 and SiO6. The structural difference of these states is only in the proportions of those blocks. At high pressure the SiO5 and SiO6 are more stable with respect to SiO4. Conversely, SiO4 is stable at low pressure. Thereby, the compression and relaxation (branches 1 and 2 in Fig. 2) result in increasing SiO5 fraction up to a maximum at 20 GPa. Further increasing Pcomp leads to convert SiO5 into SiO6. The compress– decompress procedure in present simulation differs from one reported in Ref. [22] in that we perform a long relaxation at pressure Pcomp. This relaxation allows to obtaining the high-pressure equilibrated model. Decompression and relaxation (branches 3 and 4 in

0.0 0

m

70

Pcomp (GPa) Fig. 5. The dependence of fraction of unit SiOx on pressure Pcomp.

structed models are summarized in Table 1. This result is very close to simulation data in [23] and consists with experimental data from Refs. [28,29]. As shown in Table 1 the noticeable systematic change in structure for different models concerns the height of first peak in the partial radial distribution function (PRDF) and the averaged coordination number. However, the bond length Si–O is almost unchanged and the nearest neighbor distances Si–Si, O–O slightly vary with Pcomp. This result is also clearly evidenced in Fig. 4 which demonstrates the PRDFs of considered models. The tetrahedral network structure can be seen through the coordination number Zij and the proportion of basic unit SiOx, where x = 4, 5, 6. To calculate the coordination number we use the cut-off distance, chosen as the positions of minimum after the first peak in the PRDFs. The fraction of unit SiOx is shown in Fig. 5 and one can see that the SiO4 predominates in LDA state (model M1) indicating the tetrahedral network structure. The denser model contains also an additional number of units SiO5 and SiO6 (see Fig. 5). More detail about network structure can be derived from bondangle distribution (BAD). Here we calculate only most important ones including O–Si–O and Si–O–Si angles. The first angle determines the local atomic arrangement inside unit SiOx which consti-

0.30 0.25

Pcomp= 0 GPa Pcomp=10 GPa Pcomp=20 GPa Pcomp=40 GPa Pcomp=65 GPa

(a)

0.20 0.15

Pcomp= 0 GPa Pcomp=10 GPa Pcomp=20 GPa Pcomp=40 GPa Pcomp=65 GPa

(b)

0.10

Fraction

0.05 0.00 0.25

Pcomp= 0 GPa Pcomp=10 GPa Pcomp=20 GPa Pcomp=40 GPa Pcomp=65 GPa

(c)

0.20 0.15

Pcomp= 0 GPa Pcomp=10 GPa Pcomp=20 GPa Pcomp=40 GPa Pcomp=65 GPa

(d)

0.10 0.05 0.00 30

60

90

120

150

180

60

90

120

150

180

Degree Fig. 6. The bond-angle distribution: (a) for O–Si–O in SiO4; (b) for O–Si–O in SiO5; (c) for O–Si–O in SiO6; (d) for Si–O–Si. Typical error bar is ±4.05% of data.

1219

L.T. Vinh et al. / Journal of Non-Crystalline Solids 355 (2009) 1215–1220

Fig. 2) converts partially SiO5 and almost fully SiO6 into SiO4 (see Fig. 3). The maximal density of considered models occurs at Pcomp = 20 GPa, because the glass at this pressure (branch 2) attains a maximal number of SiO5 in comparison with other compress–decompress procedure (other Pcomp); therefore, after decompressing to ambient pressure the number of remaining SiO5 is also maximal and the model is most dense. For Si–O–Si BAD we also observe two peaks: the main peak is located at 140o and a small one at 95o. As the pressure Pcomp increases, the height of main and small peaks varies in the range of 0.088–0.101 and 0.027–0.052, respectively. This data clearly indicate that the appearance of small peak in Si–O–Si BAD relates to units SiO5 and SiO6. The position of main peak changes from 130 to 147o. For model M1 that value is 145 ± 3o and it is consistent with experimental data for Si–O–Si BAD of silica glass at ambient pressure (147o [28–30]). Additional information about the linkages between two adjacent units SiOx can be inferred from the oxygen connectivity distribution. Two adjacent units SiOx are linked each to other through common oxygen atoms (see Fig. 1(f) and (e)). Table 2 shows that most of the connectivity is one or two-oxygen and their fraction strongly depend on the density of model. In particular, the fraction of two-oxygen connectivity for HDA state is 9.747% and 6 times bigger than one for LDA state.

Table 3 The VC volume distribution. The second row indicates the volume range from m VSi to (m+1) VSi. Here m = 0, 1...5 and VSi is the volume of silicon atom. The following row shows the number of VCs. For example, in model M2 there are 211 VC with the volume in the range from VSi to 2VSi. Model

Volume range

M1 M2 M3 M4 M5 M6 M7

0–1

1–2

2–3

3–4

4–5

>5

1591 1633 1859 2095 2115 1913 1878

206 211 211 151 141 171 198

60 62 34 15 14 41 33

19 18 3 2 1 9 12

6 5 0 0 0 3 4

2 1 0 0 0 0 1

3.2. Void and void aggregation Fig. 7 shows the radius void distribution (RVD) for different states of silica glass. We can see that there is a pronounced peak and about 18% of all voids in LDA model attain a radius of 0.9 ± 0.05 Å. As the density of model increases, the position of the peak shifts from 0.9 to 0.8 Å. In comparison with HDA state, the RVD shape of LDA state is wider and height of its peak is lower. To give additional insight into void aggregation, we have calculated VC volume. The VC’s volume is computed by generating 5000 points in a cube containing VC inside. Then a number of points nin located in VC is determined. The VC volume is calculated by VVC = Vcubenin/5000. Here Vcube is the volume of the cube, where 5000 points are randomly distributed. The VC volume distribution is listed in Table 3. Here we found a number of very large VCs. For example, in model M1 there are 8 VCs with volume four times bigger than 13.04 Å3 (volume of silicon atom). Obviously, that VC cannot be considered as a vacancy due to its large size in comparison with atomic size. More appropriately, it is a microscopic cavity. Such cavity is not detected in HDA state.

Fig. 8. (a) Some typical VCs; the numbers under the picture indicate the number of void in VC; (b) a part of the largest VT in a cubic box with size of 10  10  10 Å3 in the simulation cell corresponding to Pcomp of 0, 20 and 65 GPa.

60

All voids OVV SVV Largest VT

0.20 Pcomp= 0 GPa Pcomp=10 GPa Pcomp=20 GPa Pcomp=40 GPa Pcomp=65 GPa

Fraction

0.15

0.10

Volume fraction (%)

50 40 30 20 10 0

0.05

0

10

20

30

40

50

60

70

Pcomp (GPa)

0.00

Fig. 9. The dependence of volume fraction of different void kinds on pressure.

0.0

0.4

0.8

1.2

radius (Angstrom) Fig. 7. The radius void distribution.

1.6

2.0 A part of largest VT and some typical VCs are presented in Fig. 8. We also plot the volume fraction of different void kinds versus

1220

L.T. Vinh et al. / Journal of Non-Crystalline Solids 355 (2009) 1215–1220

Table 4 The characteristics of void and void aggregation. Here NV, NVSi, NVO, NVC, NVT and NVinVT are the numbers of all voids, SVVs, OVVs, VCs and VTs and the number of voids in largest VT, respectively. Model

NV

NVO

NVSi

NVC

NVT

NVinVT

M1 M2 M3 M4 M5 M6 M7

7376 7895 8069 8083 8185 8030 8056

6421 6836 6532 6035 5910 6277 6419

403 368 107 54 32 183 208

1885 1930 2107 2263 2271 2137 2126

333 367 484 613 548 504 446

5717 5808 3983 722 1353 3975 4809

Pcomp in Fig. 9. For LDA state the largest VT is spread over whole model space and for HDA state its volume significantly reduces. The great decrease in total volume of all OVVs and SVVs is also observed between HDA and LDA states. In particular, the volume fraction of SVV for model M5 is 1.2%, whereas for LDA state (model M1), this value is 15.27%, i.e. the reduction is more than 12 times. Data in Table 4 presents the number of different void kinds and it again shows that the largest VT in model M1 contains most of voids (78% all voids and 89% OVVs) and hence it is spread over whole model space. Amorphous solids may have two kinds of vacancy: thermally activated vacancy and ‘native’ vacancy [31]. First type is similar to vacancy in crystal. Their concentration depends on temperature and formation energy. The ‘native’ vacancy is associated with disordered structure and their concentration weakly depends on temperature. A great number of OVV and SVV detected in our models indicates the existence of both vacancy kinds and ensures the possibility of vacancy mechanism in silica glass. The microscopic cavity and very long VT makes diffusion mechanism in LDA system more complicated in comparison with HDA system. Silicon atom can jump into neighbor SVV and it likes the diffusion through the usual vacancy mechanism. If this movement occurs nearby microscopic cavity, then it likes the diffusion near the free surface between ‘gas’ and ‘solid’ phases. Furthermore, this movement may leads to collapse of cavity and causes the collective motion of large group of atoms. Oxygen atom has a possibility to diffuse through the long VT which is spread whole system. Probably, this kind of diffusion could take place for non-bonded oxygen and their diffusion rate is much faster than one for bonded oxygen. Such, the great difference in local structure between HDA and LDA states may lead to different diffusion mechanisms occurred in them. Now we turn the discussion on densification mechanism and the possibility of amorphous–amorphous transition. In accordance to experimental studies [12,32] the irreversible densification of SiO2 glass increases in the magnitude with increasing compress pressure up to 25 GPa, but no further densification is observed to occur for higher Pcomp. The simulation in Ref. [11] predicts the existence of HDA, LDA and intermediate states which make a partial progress from HDA to LDA state. Our simulation reproduces the experimental data in Refs. [12,32] and confirms the hypothesis suggested in [11]. Moreover, we show that the major difference in microstructure between HDA and LDA states is in different proportions of units SiO4, SiO5 and SiO6. Hence, the densification occurs with gradual variation of the proportion of those units. The

topology of units SiO4, SiO5 and SiO6 are identical for different states of silica glass. 4. Conclusions The MD simulation shows seven amorphous models SiO2 with density ranging from 2.7 to 3.15 g/cm3. These models correspond to LDA, HDA and intermediate dense states, indicating the irreversible densification of silica glass upon different compress–decompress procedure. HDA model is obtained upon compress pressure Pcomp around 20 GPa. According to our topology analysis, all amorphous states of silica consist of identical units SiO4, SiO5 and SiO6. The major difference in microstructure between HDA and LDA states is only in different relative proportion of those units. This result implies that the densification of silica glass occurs with gradual change in the proportion of units SiO4, SiO5 and SiO6. Furthermore, we find a number of cavity and very large VT in LDA model which is spread over whole system. Due to quite different void structure the diffusion mechanisms occurred in LDA and HDA states must be different. References [1] E.P. Gusev, in: G. Pacchioni, L. Skuja, D. Griscom (Eds.), Defects in SiO2 and Related Dielectrics: Science and Technology, NATO Science Series, Kluwer Academic Press, Dordrecht, 2000, p. 557 (and references therein). [2] C.R. Kurkjian, D.M. Krol, in: R.B. Devine, J.-P. Duraud, E. Dooryheé (Eds.), Structure and Imperfections in Amorphous and Crystalline Silicon Dioxide, John Wiley and Sons, New York, 2000, p. 449 (and references therein). [3] O. Mishima, L.D. Calvert, E. Whalley, Nature 310 (1984) 343. [4] O. Mishima, L.D. Calvert, E. Whalley, Nature 314 (1985) 76. [5] K.H. Smith, E. Shero, A. Chizmeshya, G.H. Wolf, J. Chem. Phys. 102 (1995) 6851. [6] Q. Williams, R.J. Hemley, M.B. Kruger, R.J. Jeanloz, Geophys. Res. 98 (1993) 22157. [7] J.-P. Iti´e, A. Polian, G. Calas, et al., Phys. Rev. Lett. 63 (1989) 398. [8] G.H. Wolf, P.F. McMillan, Structure, Dynamics and Properties of Silicate Melts, in: J.F. Stebbins, P.F. McMillan, D.B. Dingwell (Eds.), Review Mineral, vol. 32, Mineralogical Society of America, Washington, DC, 1995, p. 505. [9] A. Polian, M. Grimsditch, Phys. Rev. B 41 (1990) 6086. [10] S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, Phys. Rev. Lett. 61 (1988) 869. [11] Daniel J. Lacks, Phys. Rev. Lett. 84 (2000) 4629. [12] S. Susman, K.J. Volin, D.L. Price, et al., Phys. Rev. B 43 (1991) 1194. [13] Liping Huang, John Kieffer, Phys. Rev. B 69 (2004) 224203. [14] Liping Huang, John Kieffer, Phys. Rev. B 69 (2004) 224204. [15] Kostya Trachenko, Martin T. Dove, Phys. Rev. B 67 (2003) 064107. [16] R.L.C. Vink, G.T. Barkema, Phys. Rev. B 67 (2003) 245201. [17] James R. Rustad, David A. Yuen, Phys. Rev. A 42 (4) (1990) 2081. [18] M. Scott Shell, Pablo G. Debenedetti, et al., Phys. Rev. E 66 (2002) 011202. [19] Vo Van Hoang, Nguyen Trung Hai, Hoang Zung, Phys. Lett. A 356 (2006) 246. [20] James Badro et al., Phys. Rev. B 56 (10) (1997) 5797. [21] J.S. Tse, D.D. Klug, Y. Le Page, Phys. Rev. B 46 (1992) 5933. [22] D.J. Lacks, Phys. Rev. Lett. 84 (2000) 4629. [23] K. Vollmayr, W. Kob, K. Binder, Phys. Rev. B 54 (1996) 15808. [24] B.W.H. van Beest, G.L. Kramer, R.A. van Santen, Phys. Rev. Lett. 54 (1990) 1955. [25] P.K. Hung, H.V. Hue, L.V. Vinh, J. Non-Cryt. Sol. 352 (2006) 3332. [26] P.K. Hung, L.T. Vinh, D.M. Nghiep, P.N. Nguyen, J. Phys.: Condens. Matter 18 (2006) 9309. [27] B. Sadigh, M. Dzugutov, S.R. Elliott, Phys. Rev. B 59 (1999) 1. [28] L.I. Tatarinova, The structure of Solid Amorphous and Liquid Substances, Moscow, Nauka, 1983. [29] R.L. Mozzi, B.E. Warren, J. Appl. Cryst. 2 (1969) 164. [30] R.L. Mozzi, B.E. Warren, J. Appl. Cryst. 3 (1970) 251. [31] P.K. Hung, D.K. Belashchenko, P.N. Nguyen, Russ. Metall. 4 (1996) 155. [32] G.D. Mukherjee, S.N. Vaidya, V. Sugandhi, Phys. Rev. Lett. 87 (2001) 195501. [33] P.K. Hung, N.T. Nhan, L.T. Vinh, Modell. Simul. Mater. Sci. Eng. 17 (2009) 025003. [34] P.K. Hung, N.V. Hong, L.T. Vinh, J. Phys.: Condens. Matter 19 (2007) 466103.