# Amorphous structures and statistical sampling in first-principles molecular dynamics: The prototypical case of glassy GeSe3

Evelyne Martin, Carlo Massobrio

► To cite this version:

Evelyne Martin, Carlo Massobrio. Amorphous structures and statistical sampling in first-principles molecular dynamics: The prototypical case of glassy GeSe3. Journal of Non-Crystalline Solids, 2025, 653, pp.123415. 10.1016/j.jnoncrysol.2025.123415. hal-04944159

HAL Id: hal-04944159

https://hal.science/hal-04944159v1

Submitted on 13 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Amorphous structures and statistical sampling in first-principles molecular dynamics: The prototypical case of glassy GeSe<sub>3</sub>

Evelyne Martin<sup>●</sup>*, Carlo Massobrio

Université de Strasbourg, CNRS, Laboratoire ICube, UMR 7357, F-67037 Strasbourg, France

## ARTICLE INFO

### Keywords:

Glass structure

Molecular dynamics

Density functional theory

## ABSTRACT

The atomic structure of glassy GeSe<sub>3</sub> is obtained via first-principles molecular dynamics (FPMD) calculations by employing a simulation cell of 480 atoms. We complement and improve results previously published (Ref. Micoulaut et al., (2013)) on a smaller system (120 atoms) characterized by a marked disagreement with neutron scattering experiments on the occurrence of Ge–Ge homonuclear bonds. The present calculations show that a finite number of Ge–Ge bonds may occur on specific trajectories depending on the initial conditions selected for the quench from the liquid state. This conclusion is substantiated by the observation that the first trajectory we produced resulted in a negligible number of such homonuclear bonds, while the second did feature some of them. In terms of tetrahedral connections (either edge-sharing or corner-sharing) FPMD results appear to overestimate the edge-sharing ones, the extent of the disagreement depending on the reference experimental probe employed for the comparison (neutron scattering or NMR).

## 1. Introduction

The structure of glassy systems belonging to the Ge<sub>x</sub>Se<sub>1-x</sub> family exhibits different degrees of deviations from perfect chemical order. For a given composition, chemical order corresponds to the arrangement maximizing the number of heteronuclear bonds. In intuitive terms, structural organization in Ge<sub>x</sub>Se<sub>1-x</sub> disordered networks results from the coexistence of GeSe<sub>4</sub> units with Se–Se (x < 0.33 case) or Ge–Ge (x >0.33 case) bonds accommodating excess Se or Ge atoms not linked within tetrahedra. Along the same lines, at the stoichiometric composition only (GeSe<sub>2</sub>), one can find in principle all Ge and Se atoms within the GeSe<sub>4</sub> tetrahedron so as to avoid any miscoordinations (coordinations other than four for Ge and two for Se) and homonuclear Se–Se or Ge–Ge bonds. However, excluding a priori Ge–Ge (Se–Se) contacts for x <0.33 (x >0.33) or any of them at x = 0.33 is not adequate for binary systems characterized by a moderate difference of electronegativity between the different species, favoring the appearance of iono-covalent bonding. Stringent evidence in this direction was provided for the first time in 2000 by the technique of isotopic substitution in neutron diffraction. Measured partial pair correlation functions for glassy GeSe<sub>2</sub> revealed the presence of nearest-neighbors Ge–Se, Ge–Ge and Se–Se interatomic distances, proving the occurrence of Ge–Ge and Se–Se bonds [1]. First-principles molecular dynamics (FPMD) results concurred with the lack of perfect chemical order in glassy GeSe<sub>2</sub>, featuring various levels of agreement between

experiments and theory for the different g<sub>aβ</sub>(r) partial pair correlation functions [2–4]. Among those, g<sub>GeGe</sub>(r) was found the one with the largest standard deviation, especially for distances smaller than 4 Å, when calculated on several uncorrelated trajectories.

The absence of available experimental partial structure factors has long hampered any valuable assessment of FPMD results for Ge<sub>x</sub>Se<sub>1-x</sub> systems with composition x <0.33, as those obtained in Ref. [5] (for glassy GeSe<sub>3</sub>) and Ref. [6,7] (for glassy GeSe<sub>4</sub>). In the specific case of glassy GeSe<sub>3</sub>, results obtained for a system of 120 atoms, averaged on four independent trajectories lasting 84 ps, showed an unambiguous three peak structure for g<sup>FPMD</sup><sub>GeGe</sub>(r) corresponding to a Ge subnetwork made of Ge–Ge homonuclear bonds as well as edge-sharing and corner-sharing connections [5]. The high intensity of the peak recorded at short Ge–Ge distances points towards a strong resemblance between the pair correlation functions g<sup>FPMD</sup><sub>GeGe</sub>(r) in glassy GeSe<sub>2</sub> and GeSe<sub>3</sub>, as if the network did not undergo any major rearrangement by lowering the concentration of Ge atoms. Also, it is somewhat surprising that a deficit of Ge atoms expected to form tetrahedral connections with Se atoms did not result into an absence of homonuclear bonds when moving from GeSe<sub>2</sub> to GeSe<sub>3</sub>. In 2019 the results of Ref. [5] have been firmly questioned by a set of data based on neutron diffraction (ND) with isotopic substitution making available the partial structure factors S<sup>ND</sup><sub>GeGe</sub>(k), S<sup>ND</sup><sub>GeSe</sub>(k), S<sup>ND</sup><sub>SeSe</sub>(k) and the partial pair correlation functions g<sup>ND</sup><sub>GeGe</sub>(r), g<sup>ND</sup><sub>SeSe</sub>(r), g<sup>ND</sup><sub>GeSe</sub>(r) [8]. There is a favorable agreement between

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: evelyne.martin@cnrs.fr (E. Martin).</footnote>
these data and FPMD ones in reciprocal space, confirming the overall broad predictive power of this level of theory on short and intermediate range distances.

In this context, the following considerations are in order. For short range properties and in particular for those pertaining to chemical order (heteronuclear connections inherent in the basic structural units, such as the tetrahedron, characterized by a predominant ionic character), statistical uncertainties typical of affordable FPMD runs (say, some hundreds of ps) do not exceed a few percent, these values being lower in the presence of a more effective statistical sampling, as in liquid. Structural properties related to intermediate range order are found less accurate when compared to experiments, due to their higher sensitivity to the size of the system. When deviations from chemical order come into play, through the presence of homonuclear linkages between more mobile atoms covalently bonded, the situation can become drastically different, with larger differences between averages taken on independent trajectories. This affects the overall uncertainty of related properties and occurs exactly in the present case of glassy GeSe₃. Indeed, a remarkable similarity between Ge–Se and Se–Se partial pair correlation function is noticeable, both in terms of peak positions and intensities. However, at odds with the results of Ref. [5], $g_{GeGe}^{ND}(r)$ features no sign of homonuclear Ge–Ge bonds.

In this paper we show that, within FPMD, the simple consideration of a larger system, combined with a careful analysis of each trajectory as an entirely reliable statistical entity, allows recovering a much better agreement with experiments for a basic structural property (the pair correlation function $g_{GeGe}^{FPMD}(r)$) describing the Ge–Ge subnetwork of glassy GeSe₃. However, the very extended time scales inherent in sampling the configurational phase space in the amorphous phase remain a very hard stake to be surmounted within MD. This leads to results that depend on the initial conditions for structural details as the presence or absence of a moderate number of homonuclear bonds. In the present context, this hypothesis is confirmed when adding to a first trajectory, globally in agreement with the essence of experimental structural data, the production of a second trajectory, issued from a totally uncorrelated starting point in the liquid phase. In this latter case, the agreement with experiments is much less satisfactory in terms of occurrence of Ge–Ge homonuclear bonds, as if the trajectory had sampled a topologically different portion of the configurational landscape. Similar considerations were developed when considering the energy landscape properties of Ge–S glassy systems in connection with the rigidity transition [9].

Overall, our results drive a two-fold consideration. First, FPMD is perfectly capable to produce a Ge–Ge subnetwork of glassy GeSe₃, in agreement with ND scattering data. Second, this capability is limited by memory effects and sampling inadequacies, making the results strongly dependent on the initial conditions and/or size and cooling effects.

The paper is organized as follows. In Section 2 we provide details on our FPMD methodology. The results devoted to the structural properties (Section 3) are divided in three subsections. In the first (Section 3.1) we recall the basic definitions of partial structure factors and pair correlation functions, shown in Section 3.2 and Section 3.3, respectively. Further considerations on the network topology are developed in Section 3.4. Conclusions are drawn in Section 4.

## 2. Methodology

Our simulations were performed on a system containing $N_{at} = 480$ (120 Ge and 360 Se) atoms in a periodically repeated cubic cell of size 24.28 Å, corresponding to the experimental density of the glass at 300 K, $\rho_{exp} = 4.309$ g cm⁻³ [8].

The electronic structure was described in the framework of density functional theory (DFT) with the generalized gradient approximation (GGA) due to Becke (B) for the exchange energy and Lee, Yang and Parr (LYP) for the correlation energy [10,11]. The BLYP approach was chosen since it proved to give a better description of short-range

<figure><img src="image_3.png" /><figcaption>Fig. 1. Ionic temperature (black line) and temperature of the fictitious electronic degrees of freedom (red line) at T = 600 K.</figcaption></figure>

properties (especially in the case of Ge–Ge interactions) in a Ge–Se networks where the tetrahedral coordination is predominant (liquid and glass) [12]. In particular, we refer to Ref. [13] for several examples underlying the better performances of the BLYP approach when compared to the Perdew and Wang scheme [14]. The method by Car and Parrinello was employed to ensure a self-consistent evolution of the electronic structure during the molecular dynamics motion [15,16].

In our work, the valence electrons were treated explicitly, in conjunction with norm conserving pseudopotentials of the Trouiller–Martins type to account for core–valence interactions [17]. The wave functions were expanded at the Γ point of the supercell on a plane wave basis set with an energy cutoff $E_c = 30$ Ry. A fictitious electron mass of 1000 a.u. (i.e. in units of $m_e a_0^2$ where $m_e$ is the electron mass and $a_0$ is the Bohr radius) and a time step of $\Delta t = 0.24$ fs are adopted to integrate the equations of motion.

The strict adiabatic coupling between these two sets of degrees of freedom is ensured by controlling their temperatures via the implementation of Nosé–Hoover thermostats [18–20]. This approach proved necessary even at relatively low temperatures, due to the moderate energy separation between occupied (treated explicitly as degrees of freedom within the Car–Parrinello scheme) and unoccupied electronic states. A typical example of the adiabatic conditions achieved in our calculations is shown in Fig. 1.

A great care was exercised to prepare the glassy states issued from a starting disordered configuration featuring high diffusion followed by quenching steps down to room temperature. To achieve these conditions one has to produce a so-called thermal cycle consisting of a ramp up and a subsequent ramp down. The purpose of the ramp up is to obtain a liquid-like system at temperatures allowing effective diffusion behavior and loss of memory of the initial configuration. The ramp down aims at creating a disordered state having the feature of a glassy system (vanishing atomic movement) by quenching the liquid issued from the ramp up at affordable rates, which are typically of the order of $10^{13}$ K s⁻¹ (100 K temperature reduction followed by a trajectory of 10 ps). As a prerequisite to the implementation of the above steps, a configuration was created by exploiting a set of 480 atomic coordinates of liquid GeSe₂ produced previously [13]. Then we changed the chemical identity of 40 Ge atoms into Se atoms, so as to obtain the right concentration for the GeSe₃ composition. Optimization of the resulting atomic structure via energy minimization to T = 0 K makes available an initial starting configuration for the ramp up trajectories intended to heat the system. At this point, the heating schedule is implemented by running FPMD for 10 ps at T = 100 K, 10 ps at T = 300 K, 10 ps at T = 450 K, 20 ps at T = 600 K, 40 ps at T = 900 K, 10 ps at T = 1000 K and 20 ps at T = 1200 K. Since the total average displacement of the atomic species amounts to ~ 30 Å, with diffusion coefficients as high as $3 \cdot 10^{-5}$ cm² s⁻¹ (as extracted from the
<figure><img src="image_4.png" /><figcaption>Fig. 2. Average mean square displacements calculated for the trajectory pertaining to system A at T = 1200 K.</figcaption></figure>

asymptotic time behavior of the mean square displacements, see Fig. 2) one can safely consider that no memory of the initial configuration has been kept, the system being ready to go through the ramp down cooling trajectories. In Ref. [5] (system C hereafter), four different trajectories were created, separated by 5 ps and contributing to the statistical averages on an equal footing, as if they were fully uncorrelated.

Here we adopt a different procedure, based on a much larger interval between the starting points of the trajectories, by ensuring minimal correlation between the native configuration of the two quenched structures produced from the liquid. Accordingly, two independent sets of trajectories at different temperatures were created, with very similar lengths (system A: 30 ps at T = 1000 K, 50 ps at T = 800 K, 50 ps at T = 600 K, 30 ps at T = 450 K and 50 ps at T = 300 K; system B: 20 ps at T = 1200 K, 35 ps at T = 1000 K, 50 ps at T = 800 K, 50 ps at T = 600 K, 35 ps at T = 450 K and 50 ps at T = 300 K) and averages taken at T = 300 K on the whole interval pertaining to that temperature. The overall temperature reduction ranges from 160 to 220 ps when considering or not the final portion at T = 300 K, for a quench rate $q_r$ close to ~2 ·10¹³ K s⁻¹. This value is common to most FPMD simulations in the area of amorphous/glassy systems and it was proved to be realistic enough to avoid totally unrelaxed configurations. Both trajectories of system A and system B have been exploited to produce structural properties, with the intention to highlight their intrinsic difference that are not clearly perceptible when focusing on a global average only. This is exactly the key feature of the present work, nurtured by the observation that the intrinsic limits of FPMD calculations in terms of statistical sampling of amorphous systems can be put to good use to pinpoint the atomic structure most compatible with experimental evidence.

### 3. Structural properties

#### 3.1. Calculation procedure

Comparison with available experimental data on the structural properties of glassy GeSe₃ will be achieved by resorting to the partial pair correlation functions $g_{a\beta}(r)$ in real space and to the partial structure factors $S_{a\beta}(k)$ in reciprocal phase, expressed in the Faber-Ziman form [21,22]. These quantities are mutually related as follows:

$$g_{a\beta}(r) - 1 = \frac{1}{2\pi^2\rho r} \int_0^\infty dk \, k [S_{a\beta}(k) - 1] \sin(kr) \quad (1)$$

where $\rho$ is the atomic number density $N_{at}/V$ (V being the volume) and

$$S_{a\beta}(k) - 1 = \frac{4\pi}{k} \rho \int_0^\infty dr \, r [g_{a\beta}(r) - 1] \sin(kr). \quad (2)$$

While the isotopic substitution in neutron diffraction allows accessing each $g_{a\beta}(r)$ from the corresponding measurable $S_{a\beta}(k)$, molecular dynamics simulations target the direct calculation of $g_{a\beta}(r)$ as given by

<figure><img src="image_5.png" /><figcaption>Fig. 3. Partial structure factors of glassy GeSe₃ obtained for system A by taking the time average of Eq. (4).</figcaption></figure>

$$g_{a\beta}(r) = \frac{\langle N_{a\beta}[r \in (r_{i\alpha\beta} \pm \Delta r)] \rangle}{N_{at} 4\pi r^2 \Delta r \rho}. \quad (3)$$

In Eq. (3) $r_{i\alpha\beta}$ is the interatomic distance between a given atom i of species $\alpha$ and a given atom j of species $\beta$ and $r \in (r_{i\alpha\beta} \pm \Delta r)$ means that the values of $g_{a\beta}(r)$ are calculated via discretization of the interatomic distances in shells of width $\Delta r = 0.03-0.05$ Å. $\langle N_{a\beta} \rangle$ stands for the statistical average of the molecular dynamics trajectory. Eq. (3) is normalized so as to attain the asymptotic value of 1 for large distances, where atomic positions are totally uncorrelated.

The above expression allows for an evaluation of the partial structure factors $S_{a\beta}(k)$ via Eq. (2) devoid of the large statistical noise affecting (especially for amorphous systems) MD calculations based on the direct expression of $S_{a\beta}(k)$

$$S_{a\beta}(k) = \frac{1}{N_{at}} \sum_{i=1}^{N_{\alpha}} \sum_{j=1}^{N_{\beta}} e^{-ik \cdot (r_{i\alpha} - r_{j\beta})}. \quad (4)$$

where in Eq. (4) $N_\alpha$, $N_\beta$ are the number of atoms of species $\alpha$ and $\beta$, respectively. An example of direct calculation of $S_{a\beta}(k)$ via Eq. (4) as the average over the trajectory pertaining to system A is provided in Fig. 3. The level of noise is particularly important for $S_{GeGe}(k)$, preventing a clear-cut assessment of the level of agreement with experiments.

Over the years, the combined use of Eq. (3) (for the calculation of the partial pair correlation functions) and Eq. (2) (for the calculation of the partial structure factors) has provided a wealth of reliable results in the area of disordered and glassy systems, due to increased availability of larger periodic boxes allowing for meaningful integrations of the pair correlation functions covering the asymptotic regime at large distances.

It is worth mentioning that an alternative expression for the partial structure factors can also be found in the literature, focussing on the notion of number-number and concentration-concentration correlations, or a combination of the two. These are the so-called Bhatia-Thornton (BT) partial structure factors [23,24]

$$S_{NN}(k) = c_{\alpha} c_{\alpha} S_{aa}(k) + 2c_{\alpha} c_{\beta} S_{a\beta}(k) + c_{\beta} c_{\beta} S_{\beta\beta}(k), \quad (5)$$

$$S_{CC}(k) = c_{\alpha} c_{\beta} \{1 + c_{\alpha} c_{\beta} [(S_{aa}(k) - S_{a\beta}(k)) + (S_{\beta\beta}(k) - S_{a\beta}(k))]\}, \quad (6)$$

$$S_{NC}(k) = c_{\alpha} c_{\beta} \{c_{\alpha} (S_{aa}(k) - S_{a\beta}(k)) - c_{\beta} (S_{\beta\beta}(k) - S_{a\beta}(k))\}. \quad (7)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Faber-Ziman partial structure factors for the two trajectories considered in this work (system A and system B). The label “Exp” stands for the neutron scattering data of Ref. [8].</figcaption></figure>

The total neutron structure factor can be obtained by combining the Bhatia-Thornton (BT) partial structure factors as follows:

$$S_T(k) = S_{NN}(k) + A[S_{CC}(k)/c_\alpha c_\beta - 1] + B S_{NC}(k). \quad (8)$$

In Eq. (8), $c_\alpha$, $c_\beta$, are the atomic fractions of chemical species $\alpha$ and $\beta$, $A = c_\alpha c_\beta \Delta b^2/(\langle b \rangle)^2$, $B = 2\Delta b/(\langle b \rangle)$, $\Delta b = b_\alpha - b_\beta$, $\langle b \rangle = c_\alpha b_\alpha + c_\beta b_\beta$ is the mean coherent neutron scattering length composed with the coherent neutron scattering lengths of chemical species $\alpha$ and $\beta$, namely $b_\alpha$ and $b_\beta$.

## 3.2. Reciprocal space results: partial structure factors

Fig. 4 shows the partial structure factors $S_{GeGe}(k)$, $S_{GeSe}(k)$ and $S_{SeSe}(k)$ obtained for glassy GeSe$_3$ when calculated for system A and system B. The agreement with the experimental data of Ref. [8] is favorable in the case of $S_{GeSe}(k)$ and $S_{SeSe}(k)$, with the main features correctly reproduced at the right positions and a remarkable agreement on the intensities of the peaks. A limited discrepancy is observed for $S_{GeSe}(k)$ on the position of the first peak around ~1 Å$^{-1}$, this feature being reminiscent of the FSDP (first sharp diffraction peak), an unambiguous signature of intermediate range order. The level of agreement worsens for $S_{GeGe}(k)$, even though a broad agreement on the overall pattern of this partial structure factor does exist. Experiments and theory differ by the intensities of the peaks and some shifts in the oscillations as a function of the wavevector. This is an unmistakable sign of a reduced level of agreement to be expected also in real space, although the precise prediction of the impact of the observed difference on specific features of the partial pair correlation function remains prohibitive.

On the side of the BT partial structure factors (see Fig. 5), our result confirms the very good agreement between calculated and measured $S_{NN}(k)$, this latter quantity being a very good estimate of the total neutron structure factor $S_T(k)$ as already pointed out in Ref. [5] for system C. In fact, for networks belonging to the Ge$_x$Se$_{1-x}$ family like GeSe$_9$, GeSe$_4$, GeSe$_3$, Ge$_2$Se$_3$, GeSe$_2$, the relationship $S_T(k) \approx S_{NN}(k)$ (1% difference at most as a function of the composition) holds true

<figure><img src="image_7.png" /><figcaption>Fig. 5. Bhatia-Thornton (BT) partial structure factors for the two trajectories considered in this work (system A and system B). The label “Exp” stands for the neutron scattering data of Ref. [8].</figcaption></figure>

since $b_{Ge} = 8.185$ fm and $b_{Se} = 7.97$ fm, leading to vanishing values of A and B in Eq. (8). Interestingly, the partial structure factor $S_{CC}(k)$ does exhibit a peak in the FSDP region at ~1 Å$^{-1}$ thereby indicating the presence of fluctuations of concentration on intermediate range distance, as observed for other Ge$_x$Se$_{1-x}$ compositions [25].

## 3.3. Real space results: pair correlation functions

Based on the above analysis of the partial structure factors, both system A and system B are equally well suited to account for the atomic structure of glassy GeSe$_3$ and confirm what found for system C where the comparison with experiments did not involve the partial structure factors, at that time unavailable. However, in view of the analysis produced above for the partial structure factors, one can expect that the agreement with experiments will be higher for $g_{GeSe}(r)$ and $g_{SeSe}(r)$, in line with the consideration that the neutron scattering experiments have pointed out the absence of homonuclear Ge-Ge bonds found in system C (see Ref. [5]). In addition, it is important to underline that, due to the intrinsic character of reciprocal space data that are convolution of contributions in direct space, a conclusive insight into the network structure can only be gained by considering the pair correlation functions $g_{\alpha\beta}(r)$. These have been obtained for systems A, B and C and by neutron diffraction with isotopic substitution (Ref. [8]). The comparison is given in Fig. 6. Negligible differences among the four sets of results are recorded for $g_{GeSe}(r)$ and $g_{SeSe}(r)$, with nearest neighbor coordination number very close to those pertaining to the chemically ordered network model (CON) (see Table 1). On the contrary, the shape of $g_{GeGe}(r)$ reflects a much higher sensitivity to the model production, exhibiting a) a clear underestimate of the intensity of the third peak at ~3.7 Å (related to corner-sharing, CS, tetrahedra) common to all FPMD results b) higher intensities for the second peak at ~3 Å (representative of edge-sharing, ES, tetrahedra), especially for system B and c), most importantly, different number of homonuclear bonds at ~2.3 Å as indicated by the various heights of the first peak. The disagreement on the intensities of the second and third peak of $g_{GeGe}(r)$, common to systems A, B and C, corresponds to different proportions of ES vs CS tetrahedra. As mentioned in the introduction, more puzzling is the case of Ge-Ge contacts since experiments rule out such bonds for this system, in perfect agreement with the CON model but at odds with what obtained for system C. The present study provides clues to interpret this
<figure><img src="image_8.png" /><figcaption>Fig. 6. Partial pair correlation functions. System A, system B: FPMD trajectories produced in this work. System C: Ref. [5]. Exp: Ref. [8].</figcaption></figure>

result and bypass this apparent shortcoming. One notices that system A features a very small number of homonuclear bonds, in very close agreement with experiments. This is an unmistakable sign that the Ge–Ge subnetwork is very much sensitive to the sampling of the initial conditions, by proving that in any case no conclusions can be drawn on the basis of a single trajectory (albeit quite extended, up to 100 ps), especially for a system made of only 120 atoms, as system C. The legitimacy of this statement is boosted by the observation that system B performs much worse than system A (and somewhat also than system C, at least considering the intensity of the ES peak).

Previous work on the structure of amorphous networks (especially chalcogenides) by FPMD has pointed out that the strong variability of the Ge subnetworks in terms of nearest-neighbors arrangements can be dealt with by resorting to a large number of trajectories, allowing for the evaluation of a suitable statistical error [2]. While this holds true in any case, this strategy appears quite often unaffordable and cannot be employed systematically since the number of trajectories within the reach of FPMD is necessarily limited by computational costs. In this paper, we propose to look at the averages obtained from each individual trajectory, by the best suited to be the counterpart of the real system. This is done by acknowledging that trajectories can sample parts of the phase space more or less compatible with experimental findings.

### 3.4. Real space results: coordination numbers and tetrahedral connections

Once established the relevant pair correlation function $g_{a\beta}$ involving pairs of atoms $\alpha$ and $\beta$, one can obtain the coordination numbers $n_{\alpha}^{\beta}$ by integrating $g_{a\beta}$ up to the first minimum. The relationship to be employed to obtain the coordination numbers is the following:

$$n_{\alpha}^{\beta} = 4\pi\rho c_{\beta} \int_{r_1}^{r_2} dr r^2 [g_{a\beta}(r)] \quad (9)$$

where $r_1$ and $r_2$ are two appropriate distances expressing the range over which $n_{\alpha}^{\beta}$ is defined. By definition the total coordination numbers for atoms $\alpha$ and $\beta$ are given by $n_{\alpha} = n_{\alpha}^{\beta} + n_{\alpha}^{\beta}$ and $n_{\beta} = n_{\beta}^{\beta} + n_{\beta}^{\alpha}$ where $n_{\alpha}^{\beta}/c_{\beta} = n_{\beta}^{\alpha}/c_{\alpha}$ and the average total coordination number $n_{tot}$ is equal to $c_{\alpha} n_{\alpha} + c_{\beta} n_{\beta}$. As shown in Table 1, the three FPMD models are consistent with the notion of a chemical ordered network, with coordination numbers $n_{Ge}$, $n_{Se}$ broadly agreeing over system A, system B and system C. The most notable difference between the three sets of

<caption>**Table 1**</caption>

Values for the coordination numbers $n_{Ge}^{Ge}$, $n_{Ge}^{Se}$, $n_{Se}^{Se}$, $n_{Se}^{Ge}$, $n_{Ge}$, $n_{Se}$, and the total coordination number $n_{tot}$. Exp: neutron diffraction results of Ref. [8]. System A, B: FPMD trajectories produced in this work ($N_{tot} = 480$). System C: FPMD results of Ref. [5] ($N_{tot} = 120$). RCN stands for the random covalent network model.

<table><thead><tr><th></th><th>n<sub>Ge</sub><sup>Ge</sup></th><th>n<sub>Ge</sub><sup>Se</sup></th><th>n<sub>Ge</sub></th><th>n<sub>Ge</sub><sup>Se</sup></th><th>n<sub>Se</sub><sup>Ge</sup></th><th>n<sub>Se</sub></th><th>n<sub>tot</sub></th></tr></thead><tbody><tr><td>System A</td><td>0.02</td><td>3.96</td><td>3.98</td><td>0.68</td><td>1.32</td><td>2.00</td><td>2.50</td></tr><tr><td>System B</td><td>0.07</td><td>3.94</td><td>4.01</td><td>0.71</td><td>1.31</td><td>2.00</td><td>2.50</td></tr><tr><td>System C</td><td>0.13</td><td>3.87</td><td>4.00</td><td>0.71</td><td>1.29</td><td>2.00</td><td>2.50</td></tr><tr><td>RCN</td><td>1.6</td><td>2.4</td><td>4</td><td>1.2</td><td>0.8</td><td>2</td><td>2.5</td></tr><tr><td>CON</td><td>0</td><td>4</td><td>4</td><td>0.667</td><td>1.333</td><td>2</td><td>2.5</td></tr><tr><td>Exp.</td><td>0.0</td><td>4.00</td><td>4.00</td><td>0.70</td><td>1.333</td><td>2.00</td><td>2.50</td></tr></tbody></table>

FPMD results concern $n_{Ge}^{Ge}$, as anticipated by the intensities of the first peak in $g_{GeGe}(r)$. System A is the one best reproducing the experimental outcome featuring the absence of Ge–Ge homonuclear bonds, clearly present in both system C and to a minor extent, system B.

It is also instructive to consider the amount of Ge and Se found in edge-sharing (ES) and corner-sharing (CS) configurations (see Table 2). In the specific case of GeSe$_3$, Se atoms can also be found in Se–Se contacts, this concentration differing from the stoichiometric one by the presence of Se atoms that do not bind to Ge atoms in GeSe$_4$ tetrahedra. By definition, Ge and Se atoms in ES configurations are those found in fourfold rings, in the large majority of cases chemically ordered with alternate Ge–Se–Ge–Se bonding patterns. The three FPMD models exhibit very close amounts of Ge and Se atoms in ES configurations. However, there is an overall overestimate of this kind of connections, in line with the observation that the second peak of $g_{GeGe}(r)$ is much smaller in the experimental results of Ref. [8] when compared to the third one associated to CS configurations. In this respect, system B is the one worst performing with an unrealistic $N_{Ge}(ES)/N_{Ge}(CS)$ ratio larger than 1. This is confirmed by the shape of the Ge–Se–Ge bond angle distribution $\theta_{GeSeGe}$, exhibiting two peaks for well distinct angles corresponding to ES and CS configurations (see Fig. 7). System B is characterized by a much higher first peak, indicative of $N_{Ge}(ES)/N_{Ge}(CS) \gg 1$. On the other hand, both system A and system B predicts correctly the tetrahedral angle of 109° corresponding to the main peak of $\theta_{SeGeSe}$.

Finally, it appears that, while system C features values for $N_{Ge}(ES)/N_{Ge}$, $N_{Se}(ES)/N_{Se}$ and $N_{Ge}(ES)/N_{Ge}(CS)$ in slightly better agreement with neutron scattering experiments, system A has to be preferred in view of its negligible number of homonuclear bonds. In terms of the comparison with reported measurements, it is worth noting that the overall agreement is more favorable if one considers the data of $^{77}$Se MAS NMR experiments of Ref. [26].

### 4. Conclusions

The main motivation of this work was to reconsider the atomic structure of glassy GeSe$_3$ in light of the pieces of evidence collected by neutron diffraction with isotopic substitution (Ref. [8]) making available detailed information on the atomic Ge and Se subnetworks. Previous FPMD results have provided evidence on the occurrence of Ge–Ge bonds, in striking contrast with the experiments quoted above. More specifically, Ge–Ge homonuclear bonds have appeared on the Se-rich side of the $Ge_{x}Se_{1-x}$ family of glasses in FPMD studies, in disagreement with their absence clearly demonstrated for GeSe$_3$ and GeSe$_4$ in Ref. [8]. Given the above context, we intended to ascertain whether this feature is an unavoidable outcome of FPMD approaches or it can be differently accounted for. To this goal, we produced two fully uncorrelated trajectories and quenched them from the liquid state at rates typical of FPMD calculations. In the first case, we obtained an overall very good agreement with neutron diffraction data via the analysis of the Ge–Ge pair correlation function, the number of Ge–homonuclear bonds being essentially negligible. In the second case, this same quantity is by far not negligible, pointing out a drastic change with respect to the behavior of the first trajectory, the agreement on the
<caption>Table 2</caption>

$N_{Ge}(ES)$, $N_{Se}(ES)$, $N_{Ge}(CS)$, $N_{Se}(CS)$ are the number of Ge or Se atoms in edge-sharing and corner-sharing configurations, respectively. $N_{Ge}$, $N_{Se}$ are the total number of Ge, Se atoms in system A, B or C. Exp1: Neutron diffraction results of Ref. [8]. Exp2: $^{77}$Se MAS NMR experiments of Ref. [26] as reported in Ref. [8]. Note that the value of $N_{Ge}(ES)/N_{Ge}(CS)$ given in Ref. [5] (0.32) is erroneous since not consistent with the number of Ge atoms found in the ES configuration (Fig. 11 of Ref. [5]) and with the value of $n_{Ge}^{Ge}$ (Table II of Ref. [5]).

<table><thead><tr><th></th><th>System A</th><th>System B</th><th>System C</th><th>Exp1</th><th>Exp2</th></tr></thead><tbody><tr><td>$N_{Ge}(ES)$</td><td>39 ($N_{Ge}$ = 120)</td><td>59 ($N_{Ge}$ = 120)</td><td>8 ($N_{Ge}$ = 30)</td><td></td><td></td></tr><tr><td>$N_{Ge}(ES)/N_{Ge}$</td><td>0.32</td><td>0.49</td><td>0.27</td><td>0.15</td><td>0.30</td></tr><tr><td>$N_{Ge}(CS)$</td><td>80 ($N_{Ge}$ = 120)</td><td>53 ($N_{Ge}$ = 120)</td><td>18 ($N_{Ge}$ = 30)</td><td></td><td></td></tr><tr><td>$N_{Se}(ES)$</td><td>40 ($N_{Se}$ = 360)</td><td>60 ($N_{Se}$ = 360)</td><td>8 ($N_{Se}$ = 90)</td><td></td><td></td></tr><tr><td>$N_{Se}(ES)/N_{Se}$</td><td>0.11</td><td>0.17</td><td>0.09</td><td>0.05</td><td>0.10</td></tr><tr><td>$N_{Ge}(ES)/N_{Ge}(CS)$</td><td>0.49</td><td>1.11</td><td>0.44</td><td>0.18</td><td>0.43</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Bond angle distributions $\theta_{SeGeSe}$ and $\theta_{GeSeGe}$ for system A and system B.</figcaption></figure>

shape of the other pair correlation functions remaining very favorable.

Since we found a large uncertainty on a very specific quantity, a distinction in this regard has to be made among properties pertaining to short range chemical order and those describing deviations from chemical order, such as the existence of Ge-Ge connections discussed above. While the first are found to agree with experimental data within a few percents, the latter might have much larger statistical errors, due to their dependence on the initial configurations that tend to persist during the quench procedure. Chemical order is the key feature to be considered in this context, since deviations from these conditions can involve bond breaking and restoring processes hard to describe within the time span of FPMD calculations. In view of these thoughts, the observed dependence of the number of Ge-Ge connection on the trajectory employed to produce statistical averages can be ascribed to drastic differences in the starting liquid configurations. By comparison, the size of the simulation cell has a smaller effect on this outcome, as proved by the fact that one of the trajectories we created (system B, $N_{at}$= 480) has a number of Ge-Ge connections not too dissimilar from the one of system C, for $N_{at}$= 120. Incidentally, it has to be observed that the size of the system issue is more troublesome when the number of atoms for a given species is intrinsically small, thereby affecting the overall statistical sampling.

More than describing this behavior by attributing an error bar to the mean result obtained from the two independent trajectories, we decided to consider separately the outcome of each single trajectory. This allows exemplifying the strong dependence on the initial conditions of the FPMD atomic glassy structures of GeSe$_3$. In this way we showed that the specific discrepancy between FPMD results and neutron diffraction regarding the impact of homonuclear Ge-Ge bonds data cannot be ascribed to some shortcoming of the FPMD theoretical framework.

Instead, it reflects an intrinsic limit of statistical sampling that is very hard to surmount when producing disordered systems configurationally frozen from the liquid. Having established that FPMD can very well describe glassy GeSe$_3$ without predicting inevitably homonuclear bonds, we did observe some persisting disagreement between theory and experiments on the relative percentage of edge-sharing vs corner-sharing configurations. Edge-sharing configurations are more present in FPMD models than in neutron scattering experiments while the comparison improves when considering NMR results.

We reiterate that the strategy consisting in the production of several trajectories with average values assigned to the given quantities and related error bar extracted from the standard deviation of the results is perfectly legitimate and can also be employed to describe the structural properties of glassy systems. However, it has two main disadvantages, these being the high computational cost and the unavoidable limited sampling of the configurational phase space hampering a clear distinction between behaviors differently consistent with the experimental outcome. In the present case, the average result for the homonuclear Ge-Ge bonds calculated on trajectories A and B turns out to be 0.045 with an error bar equal to 0.025. This results by itself is much less instructive than the observation and careful analysis of distinct trajectories carried out in this work.

In conclusion, a great care should be exercised when commenting on supposedly remarkable differences between theory (FPMD calculations) and experiments targeting the atomic structure of amorphous systems. The various disagreements found on various quantities can be identified as belonging to two categories: to the first pertain those due to limited statistical sampling while the second features differences more profoundly related to the description of bonding and, possibly, the impact of size effects. This calls for the use of modern data mining techniques aimed at the construction of reliable interatomic potentials opening the way to the consideration of much larger systems. Efforts in this direction are currently in progress along the lines pioneered in a recent paper [27].

### CRediT authorship contribution statement

**Evelyne Martin:** Investigation, Formal analysis. **Carlo Massobrio:** Writing – original draft, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

We acknowledge GENCI (Grand Equipement National de Calcul Intensif), France (Grant No. A0xx0905071, A0xx0910296) and the High Performance Computing Center of the University of Strasbourg, France for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data. Instructive discussions with Mauro Boero and Guido Ori are gratefully acknowledged.
## Data availability

The raw/processed data required to reproduce these findings are available from the corresponding author on reasonable request.

## References

[1] Ingrid Petri, Philip S. Salmon, Henry E. Fischer, Defects in a disordered world: The structure of glassy GeSe₂, Phys. Rev. Lett. 84 (2000) 2413-2416, http://dx.doi.org/10.1103/PhysRevLett.84.2413, URL https://link.aps.org/doi/10.1103/PhysRevLett.84.2413.

[2] Carlo Massobrio, Alfredo Pasquarello, Short and intermediate range order in amorphous GeSe₂, Phys. Rev. B 77 (2008) 144207, http://dx.doi.org/10.1103/PhysRevB.77.144207, URL https://link.aps.org/doi/10.1103/PhysRevB.77.144207.

[3] Mark Cobb, D.A. Drabold, R.L. Cappelletti, Ab initio molecular-dynamics study of the structural, vibrational, and electronic properties of glassy GeSe₂, Phys. Rev. B 54 (1996) 12162-12171, http://dx.doi.org/10.1103/PhysRevB.54.12162, URL https://link.aps.org/doi/10.1103/PhysRevB.54.12162.

[4] De Nyago Tafen, D.A. Drabold, Realistic models of binary glasses from models of tetrahedral amorphous semiconductors, Phys. Rev. B 68 (2003) 165208, http://dx.doi.org/10.1103/PhysRevB.68.165208, URL https://link.aps.org/doi/10.1103/PhysRevB.68.165208.

[5] M. Micoulaut, A. Kachmar, M. Bauchy, S. Le Roux, C. Massobrio, M. Boero, Structure, topology, rings, and vibrational and electronic properties of Ge<sub>x</sub>Se<sub>1-x</sub> glasses across the rigidity transition: A numerical study, Phys. Rev. B 88 (2013) 054203, http://dx.doi.org/10.1103/PhysRevB.88.054203, URL https://link.aps.org/doi/10.1103/PhysRevB.88.054203.

[6] Kateryna Sykina, Eric Furet, Bruno Bureau, Sébastien Le Roux, Carlo Massobrio, Network connectivity and extended se chains in the atomic structure of glassy GeSe<sub>4</sub>, Chem. Phys. Lett. 547 (2012) 30-34, http://dx.doi.org/10.1016/j.cplett.2012.07.077, URL https://www.sciencedirect.com/science/article/pii/S0009261412008986.

[7] Carlo Massobrio, Massimo Celino, Philip S. Salmon, Richard A. Martin, Matthieu Micoulaut, Alfredo Pasquarello, Atomic structure of the two intermediate phase glasses SiSe<sub>4</sub> and GeSe<sub>4</sub>, Phys. Rev. B 79 (2009) 174201, http://dx.doi.org/10.1103/PhysRevB.79.174201, URL https://link.aps.org/doi/10.1103/PhysRevB.79.174201.

[8] Ruth F. Rowlands, Anita Zeidler, Henry E. Fischer, Philip S. Salmon, Structure of the intermediate phase glasses GeSe<sub>3</sub> and GeSe<sub>4</sub>: The deployment of neutron diffraction with isotope substitution, Front. Mater. 6 (2019) 133, http://dx.doi.org/10.3389/fmats.2019.00133, URL https://www.frontiersin.org/article/10.3389/fmats.2019.00133.

[9] S. Chakraborty, P. Boolchand, M. Micoulaut, Structural properties of Ge-S amorphous networks in relationship with rigidity transitions: An ab initio molecular dynamics study, Phys. Rev. B 96 (2017) 094205, http://dx.doi.org/10.1103/PhysRevB.96.094205, URL https://link.aps.org/doi/10.1103/PhysRevB.96.094205.

[10] A.D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38 (1988) 3098-3100, http://dx.doi.org/10.1103/PhysRevA.38.3098, URL https://link.aps.org/doi/10.1103/PhysRevA.38.3098.

[11] Chengteh Lee, Weitao Yang, Robert G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37 (1988) 785-789, http://dx.doi.org/10.1103/PhysRevB.37.785, URL https://link.aps.org/doi/10.1103/PhysRevB.37.785.

[12] Matthieu Micoulaut, Rodolphe Villeumier, Carlo Massobrio, Improved modeling of liquid GeSe₂: Impact of the exchange-correlation functional, Phys. Rev. B 79 (2009) 214205, http://dx.doi.org/10.1103/PhysRevB.79.214205, URL https://link.aps.org/doi/10.1103/PhysRevB.79.214205.

[13] Carlo Massobrio, The Structure of Amorphous Materials using Molecular Dynamics, IOP Publishing, 2022, http://dx.doi.org/10.1088/978-0-7503-2436-6.

[14] John P. Perdew, Yue Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Phys. Rev. B 45 (1992) 13244-13249, http://dx.doi.org/10.1103/PhysRevB.45.13244, URL https://link.aps.org/doi/10.1103/PhysRevB.45.13244.

[15] R. Car, M. Parrinello, Unified approach for molecular dynamics and density-functional theory, Phys. Rev. Lett. 55 (1985) 2471-2474, http://dx.doi.org/10.1103/PhysRevLett.55.2471, URL https://link.aps.org/doi/10.1103/PhysRevLett.55.2471.

[16] Dominik Marx, Jürg Hutter, Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods, Cambridge University Press, 2009, http://dx.doi.org/10.1017/CBO9780511609633.

[17] N. Troullier, José Luís Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43 (1991) 1993-2006, http://dx.doi.org/10.1103/PhysRevB.43.1993, URL https://link.aps.org/doi/10.1103/PhysRevB.43.1993.

[18] Shuichi Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1) (1984) 511-519, URL https://doi.org/10.1063/1.447334.

[19] William G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31 (1985) 1695-1697, http://dx.doi.org/10.1103/PhysRevA.31.1695, URL https://link.aps.org/doi/10.1103/PhysRevA.31.1695.

[20] Peter E. Blöchl, M. Parrinello, Adiabaticity in first-principles molecular dynamics, Phys. Rev. B 45 (1992) 9413-9416, http://dx.doi.org/10.1103/PhysRevB.45.9413, URL https://link.aps.org/doi/10.1103/PhysRevB.45.9413.

[21] T.E. Faber, J.M. Ziman, A theory of the electrical properties of liquid metals, Philos. Mag. A 11 (109) (1965) 153-173, URL https://doi.org/10.1080/14786436508211931.

[22] Y. Waseda, The Structure of Non-Crystalline Materials — Liquids and Amorphous Solids, McGraw-Hill, New York, 1980.

[23] A.B. Bhatia, D.E. Thornton, Structural aspects of the electrical resistivity of binary alloys, Phys. Rev. B 2 (1970) 3004-3012, http://dx.doi.org/10.1103/PhysRevB.2.3004, URL https://link.aps.org/doi/10.1103/PhysRevB.2.3004.

[24] Philip S. Salmon, The structure of molten and glassy 2:1 binary systems: an approach using the Bhatia-Thornton formalism, Proc. R. Soc. Lond. A 437 (1992) 591, http://dx.doi.org/10.1103/PhysRevLett.103.157801, URL https://royalsocietypublishing.org/doi/10.1098/rspa.1992.0081.

[25] Ian T. Penfold, Philip S. Salmon, Structure of covalently bonded glass-forming melts: A full partial-structure-factor analysis of liquid GeSe<sub>3</sub>, Phys. Rev. Lett. 67 (1991) 97-100, http://dx.doi.org/10.1103/PhysRevLett.67.97, URL https://link.aps.org/doi/10.1103/PhysRevLett.67.97.

[26] E.L. Gjersing, S. Sen, B.G. Aitken, Structure, connectivity, and configurational entropy of Ge<sub>x</sub>Se<sub>1-x</sub>, glasses: Results from <sup>77</sup>Se MAS NMR spectroscopy, J. Phys. Chem. C 114 (18) (2010) 8601-8608, URL https://doi.org/10.1021/jp1014143.

[27] Pankaj Rajak, Nitish Baradwaj, Ken-ichi Nomura, Aravind Krishnamoorthy, Jose P. Rino, Kohei Shimamura, Shogo Fukushima, Fuyuki Shimojo, Rajiv Kalia, Aiichiro Nakano, Priya Vashishta, Neural network quantum molecular dynamics, intermediate range order in gese<sub>2</sub>, and neutron scattering experiments, J. Phys. Chem. Lett. 12 (25) (2021) 6020-6028, URL https://doi.org/10.1021/acs.jpclett.1c01272.