Xiliang Lian
email: xiliang.lian@sorbonne-universite.fr

Mathieu Salanne
email: mathieu.salanne@sorbonne-universite.fr

Capturing the interactions in the BaSnF 4 ionic conductor: comparison between a machine-learning potential and a polarizable force field

BaSnF 4 is a prospective solid state electrolyte for fluoride ion batteries. However, the diffusion mechanism of the fluoride ions remains difficult to study, both in experiments and in simulations. In principle, ab initio molecular dynamics could allow to fill this gap, but this method remains very costly from the computational point of view. Using machine learning potentials is a promising method that can potentially address the the accuracy issues of classical empirical potentials while maintaining high efficiency.

In this work, we fitted a dipole polarizable ion model and trained machine learning potential for BaSnF 4 and made comprehensive comparisons on the ease of training, accuracy and efficiency. We also compared the results with the case of a simpler ionic system (NaF). We show that contrarily to the latter, for BaSnF 4 the machine learning potential offers much higher versatility. The current work lays foundations for the investigation of fluoride ion mobility in BaSnF 4 and provides insight on the choice of methods for atomistic simulations.

I. INTRODUCTION

Fluoride ion batteries (FIBs) are a promising alternative to Lithium-ion batteries due to their improved safety, low cost, and high theoretical volumetric energy density. They have therefore become increasingly investigated in recent years [START_REF] Xiao | The case for fluoride-ion batteries[END_REF][START_REF] Nowroozi | Fluoride ion batteries -past, present, and future[END_REF][START_REF] Liu | Rechargeable anion-shuttle batteries for low-cost energy storage[END_REF] . However, the low ionic mobility of fluoride ions within the solid electrolyte and the requirement of high-temperature operation are the main bottlenecks to overcome before a potential commercialization of FIBs [START_REF] Xiao | The case for fluoride-ion batteries[END_REF] . The best performance in terms of ionic conductivity are obtained with PbSnF 4 that can reach as high as 1 × 10 -3 S cm -1 at room temperature but its application is hampered by the fact that both Pb and Sn can be oxidized and Pb will pose environmental issues.

Substituting Pb with Ba results in BaSnF 4 which exhibits reasonably high room temperature ionic conductivity, yet it remains far below that of Li-ion ionic conductors. The tetragonal phase of BaSnF 4 possesses ionic conductivity almost two orders of magnitude higher than the cubic one [START_REF] Ahmad | Structure, ionic conduction, and giant dielectric properties of mechanochemically synthesized BaSnF 4[END_REF] , and is of particular interest for solid-state electrolytes. The structure of tetragonal phase BaSnF 4 (simply referred as BaSnF 4 in the remaining of the text) is shown on Fig. 1 (a) and can be described as alternating Ba-Ba and Sn-Sn double layers. Notably, the Sn atoms have stereoactive lone pairs pointing towards the near-empty Sn-Sn layer [START_REF] Dénès | BaSnF 4 -A new fluoride ionic conductor with the α-PbSnF 4 structure[END_REF] .

Previous work using NMR technique has shown that two different fluoride ion sublattices exist [START_REF] Chaudhuri | Resolving the different dynamics of the fluorine sublattices in the anionic conductor BaSnF 4 by using high-resolution MAS NMR techniques[END_REF] . The fluoride ions lying between the Ba-Ba layers are immobile and require high energy to migrate while the fluoride ions located between the Ba-Sn / Sn-Sn (the latter with fractional occupation, not shown in the figure) layers are highly mobile and are under rapid exchange even at room temperature. Interestingly, at high temperatures, the participation of fluoride ions between Ba-Ba layers to the conduction can be observed from the 2D NMR exchange spectra. Despite a few experimental and theoretical works on BaSnF 4 , the ionic conduction mechanism between Ba-Sn/Sn-Sn layers is still unclear and understanding it would significantly benefit to the rational design of FIB electrolytes.

Modeling accurately the interactions in BaSnF 4 using computer simulations is challenging since the lone-pair interactions are tricky to capture and the relatively low ionic conductivity requires long molecular dynamics (MD) simulations. Density functional theory (DFT)-based ab initio molecular dynamics (AIMD) provides an accurate description of the interactions, but this method is limited to small system sizes and to trajectories of hundreds of picoseconds due to its high computational cost. Classical MD allows to obtain trajectories over much 2 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. longer time scales and has shown exceptional performance in numerous previous researches of ionic conductors. This is because it involves classical interaction potentials (IP) which are generally based on simple analytical expressions that try to mimic the physical phenomenon at play. In particular, previous works have proved the accuracy of polarizable models (the so-called dipole-polarizable ionic model (DIPPIM)) for simulating ionic materials, [START_REF] Madden | From first-principles to material properties[END_REF][START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF] and even ionic conductors such as PbF 2 [START_REF] Castiglione | Fluoride ion disorder and clustering in superionic PbF 2[END_REF] or Li 7 La 3 Zr 2 O 12 [START_REF] Burbano | Sparse cyclic excitations explain the low ionic conductivity of stoichiometric Li 7 La 3 Zr 2 O 12[END_REF] . However, the existence of Sn lone pairs is hard to account for when using such classical IP.

PLEASE CITE THIS ARTICLE AS

This limitation may be overcome by using data-driven approaches instead. In particular, machine learning potentials (MLPs), since the initial work by Behler and Parrinello, have received significant amount of attention [11][START_REF] Behler | Perspective: Machine learning potentials for atomistic simulations[END_REF][START_REF] Unke | Machine Learning Force Fields[END_REF] . They are expected to exhibit ab initio accuracy at a computational efficiency comparable to that of classical IPs. MLP describes the potential energy surface by establishing a statistical relationship between the energy and structures, which are described by descriptors in MLPs. These descriptors are constructed in such a way that the energy of a structure is invariant in translation, rotation and permutation of two chemically identical atoms. Previous work has shown that MLPs works well for simple systems such as Si, Ni, Cu etc. [START_REF] Zuo | Performance and Cost Assessment of Machine Learning Interatomic Potentials[END_REF] and also provide accurate results

for relatively complex materials such as Li 7 La 3 Zr 2 O 12 [START_REF] Kim | Flexible machine-learning interatomic potential for simulating structural disordering behavior of Li 7 La 3 Zr 2 O 12 solid electrolytes[END_REF] and liquids such as molten NaCl [START_REF] Li | Development of robust neural-network interatomic potential for molten salt[END_REF] .

The almost perfect agreement between MLP and DFT results reveals that MLP can capture the local environmental variations very well although the physical terms are not explicitly included. Although MLP has demonstrated excellent performance for CaF 2 17 , there has not been significant work on MLP for fluoride ion conductors and it has so far not been applied

to systems with delicate interactions as BaSnF 4 .

In this work, we fitted classical IP and trained MLP for BaSnF 4 . To make effective comparisons, we have included NaF as a reference system since it adopts a much simpler rocksalt structure (Fig. 1 (b)) for which the DIPPIM is expected to perform well. For both systems, the two methods are compared in terms of parameterization difficulty (i.e.

the amount of necessary data in the parameterization process), computational cost, and accuracy for reproducing structural, thermodynamic and vibrational properties. We show that the versatility of MLPs makes them more adequate for simulating systems with complex interactions, despite a lower computational efficiency. within the density-functional theory (DFT) framework [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]19 . The electronic wave functions were expanded in a plane-wave basis set. The core-valence interactions are described by projector-augmented wave (PAW) pseudopotentials. For all the AIMD calculations, the Brillouin zone is sampled at gamma point only to increase the efficiency and the electronic convergence threshold is 10 -6 eV. The Verlet algorithm is used to integrate Newton's equations of motion, at a time step of 1 fs on an NVT ensemble where we fix the number of atoms, the volume, and the temperature. Nosé-Hoover thermostat is used to control the temperature at 298 K. The same set of parameters were used when applying MLP.

For NaF, the exchange-correlation energy is approached with Perdew-Burke-Ernzerhof (PBE) functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] with a kinetic energy cutoff of 700 eV. The valence configuration for Na 4

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. and F are 2s 2 2p 6 3s 1 and 2s 2 2p 5 . We took the initial configuration from materials project database [START_REF] Jain | Commentary: The materials project: A materials genome approach to accelerating materials innovation[END_REF] and the AIMD calculations are made with a 2×2×2 supercell. For BaSnF 4 , the calculations are performed with PBE for solids (PBESol) exchange-correlation functional [START_REF] Perdew | On-the-fly machine learning force field generation: Application to melting points[END_REF] , which was found to best reproduce the experimental data as shown in the Supplementary Table I. The plane wave cutoff of BaSnF 4 was set to 520 eV. The valence shell configurations for Ba, Sn, and F are 5s 2 5p 6 6s 2 , 4d 10 5s 2 5p 2 , and 2s 2 2p 5 respectively. The initial configuration is taken from the experiment and a 3×3×1 supercell is constructed for AIMD calculation.

PLEASE CITE THIS ARTICLE AS

B. Machine Learning Potential

Model

We employ on-the-fly MLPs following the methodology described and implemented in VASP to generate MLPs for BaSnF 4 and NaF. The underlying theory and training parameters are extensively described in the previous works 23, [START_REF] Jinnouchi | Phase Transitions of Hybrid Perovskites Simulated by Machine-Learning Force Fields Trained on the Fly with Bayesian Inference[END_REF] . In this context, we provide a succinct overview of the pertinent theoretical aspects employed in this study only.

The potential energy surface of a given structure with N atoms is decomposed into contributions from individual atomic energies:

U = N i=1 U i
The local energy is determined by the local environments, which can be represented by the local atomic distribution:

ρ i (r) = N j=1 f cut (r ij )g(r -r ij )
where f cut is a cutoff function which excludes the contributions of environments outside the cutoff radius and g(r -r ij ) is a Gaussian-type delta function. In VASP implementation, radial (ρ

i ) and angular descriptors (ρ

i ) that can satisfy the rotational invariance are constructed to describe the local environments. Therefore, we have:

U i = F [ρ (2) i , ρ (3) i ]
The radial descriptor describes the probability of finding another atom j from the ith atom with a distance r. The angular descriptor in this context refers to the probability of 5 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. observing atom j at a distance r from atom i, and simultaneously observing atom k at a distance s from atom i, with an angle ̸ kij = θ between them. In practice, the descriptors are discretized and expanded into terms of radial basis functions and spherical harmonics:

PLEASE CITE THIS ARTICLE AS

ρ (2) i (r) = 1 √ 4π N 0 R n=1 c i n χ n0 (r) ρ (3) i (r, s, θ) = Lmax l=0 N l R n=1 N l R µ=1 2l + 1 2 p i nµl χ nl (r)χ µl (s)P l (cosθ)
and the expansion coefficients for angular descriptors are given in P i nlµ :

P i nµl = 8π 2 2l + 1 l m=-l [c i nlm c i * µlm - Na j c ij nlm c ij * µlm ]
where χ µl and P l denote normalized spherical Bessel functions and Legendre polynomials of order l respectively. Here, c i nlm , N 0 R , and N l R respectively represent the expansion coefficients, number of radial basis functions for radial descriptors ρ

x (2) T i = (c i 1 , c i 2 , ...) x (3) T i = (p i 110 , p i 111 , ..., p i 120 , p i 121 , ...)
Finally, the energy can be written as weighted kernel function K(x i , x iB ) where x i and

x iB respectively represent the local instantaneous and the local reference configurations:

F (x i ) = N B i B =1 ω iB K(x i , x iB )
The kernel function takes the following form:

K(x i , x iB ) = [β (2) (x (2) i • x (2) iB ) + β (3) (x (3) i • x(3) iB ))] ζ (3)
Here, the x(3) i denotes a normalized vector of x

i . β [START_REF] Nowroozi | Fluoride ion batteries -past, present, and future[END_REF] and β [START_REF] Nowroozi | Fluoride ion batteries -past, present, and future[END_REF] are weighting parameters. ζ [START_REF] Liu | Rechargeable anion-shuttle batteries for low-cost energy storage[END_REF] controls the sharpness of the kernel function.

6

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

Training

In this study, the MLP training process was performed on-the-fly, which involved constructing the potential during the AIMD calculations. The parameters for controlling the AIMD part of the calculations were set up as described in the previous section. Now, we provide a brief summary of the model training procedures and preliminary results.

To generate the MLP datasets for NaF, a 2×2×2 supercell containing 64 atoms was utilized. The initial structure was equilibrated at 300 K for 30 ps using AIMD before the training. Then, The ML model was trained by gradually heating the system from 250 K to 600 K, which allows for a broader sampling of the potential energy surface. In the case of BaSnF 4 , a slightly different setup was employed. A 3×3×1 supercell consisting of 108 atoms was used, and the model was trained by gradually heating the system from 200 K to 500 K over 100 ps. The initial equilibration was performed at 298 K for 50 ps using the on-the-fly training scheme. Following this equilibration, the MLP was abandoned and the final configuration obtained from the equilibration step served as the starting point for further training.

Here, we discuss some important hyperparameters used across the training. The cutoff radius for radial ρ

C. Classical Interaction Potential: DIPPIM

Model

The DIPPIM has been previously reported and successfully applied to the study of many ionic materials, including solid-state Li-ion batteries electrolytes [START_REF] Burbano | Sparse cyclic excitations explain the low ionic conductivity of stoichiometric Li 7 La 3 Zr 2 O 12[END_REF] and fluoride ionic conductors. While interaction potentials were previously reported for NaF [START_REF] Salanne | Heat-transport properties of molten fluorides: Determination from first-principles[END_REF] , they were generally aimed at simulating the molten state, so that in the current work we have fitted a new set of parameters to ensure its consistency with BaSnF 4 , a system for which the interaction potential has never been previously reported. The DIPPIM comprises four components, including charge-charge interaction, short-range repulsion, dispersion interaction, and the dipolar interaction.

The coulombic or charge-charge interactions are defined by

V qq (r ij ) = i≤j q i q j r ij
where r ij denotes the interatomic distances and q i and q j are charges carried by atoms i and j. The potential includes a short-range repulsive term with an exponential behavior within the interionic separation range. A Gaussian function is added as a steep repulsive wall to avoid unstabilities at high temperatures. This term is especially beneficial for systems with highly polarizable ions, like oxygen, as it prevents numerical instabilities in dipoles at very small anion-cation separations. Here, we have imposed the repulsive wall on Ba-F and Sn-F interactions. This short-range potential reads:

V rep (r ij ) = i≤j B ij exp(-η ij r ij ) r n ij + i≤j B ′ ij exp(-η ′ ij r 2 ij )
The dispersion interactions include dipole-dipole and dipole-quadruple terms:

V disp (r ij ) = - i≤j f ij 6 (r ij ) C 6 ij r 6 ij + f ij 8 (r ij ) D 8 ij r 8 ij
where C ij 6 and D ij 8 are the dipole-dipole and dipole-quadruple dispersion coefficients. The Tang-Toennies functions [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF] are used to damp the dispersion terms at short range: 8 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

f ij n (r ij ) = 1 -exp(-b ij n r ij ) n k=0 (b ij n r ij ) k k!
Finally, the polarization part of the potential incorporates dipolar effects:

V pol (r ij ) = i,j (q i µ j,α f ij 4 (r ij ) -q j µ i,α f ji 4 (r ij ))T (1) α (r ij ) - i,j µ i,α µ j,β T (2)
αβ (r ij ) + i 1 2α i |µ i | 2
In the above expression, α i represents the polarizability of ion i, while µ i denotes the dipoles. Additionally, T [START_REF] Xiao | The case for fluoride-ion batteries[END_REF] and T [START_REF] Nowroozi | Fluoride ion batteries -past, present, and future[END_REF] correspond to the charge-dipole and dipole-dipole interaction tensors, respectively. The components of these tensors are given by:

T (1) α (r ij ) = -r α /r 3 T (2) αβ (r ij ) = (3r α r β -r 2 δ αβ )/r 5
The influence of short-range induction effects on the dipoles is considered by utilizing the Tang-Toennies damping functions (f ij 4 ) [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF] . These damping functions, similar to those employed for damping dispersion interactions, take into account the overlap of charge densities that affect the induced dipoles. The parameter b ij n determines the range over which this overlap has an impact. Furthermore, the damping of these induction effects necessitates the inclusion of an additional pre-exponential parameter, c ij , which determines the strength of the ion's response to this phenomenon:

f ij n (r ij ) = 1 -c ij exp(-b ij n r ij ) n k=0 (b ij n r ij ) k k!
Note that in the DIPPIM, the induced dipoles are treated as a set of additional degrees of freedom, which are obtained by minimizing the energy at each timestep. This is formally equivalent to determining them self-consistently, so that the polarization interactions have a many-body character.

Fitting of the parameters and simulations details

The DIPPIM for NaF and BaSnF 4 was fitted using a generalized force-fitting method.

Initially, 10 configurations were randomly selected from the AIMD trajectory, for which 9

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. DFT calculations were performed to compute the forces. In addition, Wannier90 code was utilized to identify the Wannier centers and calculate the dipoles on each atom [START_REF] Pizzi | Wannier90 as a community code: New features and applications[END_REF] . Fitting those in addition to the individual forces allows to reduce the risk of error cancellations between the various parts of the potential [START_REF] Aguado | Multipoles and interaction potentials in ionic materials from planewave-DFT calculations[END_REF] . It is important to note that in the DIPPIM, the atomic charges are not treated as adjustable parameters; they are taken as the nominal charges of the ions, i.e. +1 for Na + , -1 for F -and +2 for Ba 2+ and Sn 2+ . The fitting process is performed by minimizing the following objective function:

PLEASE CITE THIS ARTICLE AS

χ 2 = N i |A DF T i -A IP i | 2 N i |A DF T i | 2
Here, A DF T i and A IP i comprise of the forces and dipoles for all atoms obtained from DFT calculations and the fitted interaction potential, respectively. For NaF and BaSnF 4 , the resulting objective function values were found to be 0.003 and 0.23, respectively. In the case of BaSnF 4 , the objective function could not be further reduced primarily due to the stability of the interaction potential.

All our calculations made with DIPPIM have been performed using Metalwalls [START_REF] Marin-Laflèche | MetalWalls: A classical molecular dynamics software dedicated to the simulation of electrochemical systems[END_REF][START_REF] Coretti | MetalWalls: Simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF] . The room temperature is maintained by a Nosé-Hoover chain thermostat with a characteristic period of 1 ps. Apart from this, the parameters are set to the same as those in AIMD calculations to ensure consistency.

III. RESULTS AND DISCUSSION

A. Validation of the potentials on DFT forces

To validate the two different methods on a similar dataset, we randomly select 200 structures from an independent molecular dynamics trajectory, ensuring that these configurations are proportional to the number of data points. In the case of MLP, the histogram displays a highly Gaussian distribution of the errors, showing that the errors are random. For DIPPIM, some deviations from Gaussian behavior are observed for the atoms with large absolute value of the DFT forces, which reflects the fact that the model does not fully catch the electronic structure effects. In terms of absolute error, the comparison of the root mean squared error (RMSE) reveals that the MLP slightly outperforms the DIPPIM.

In the case of BaSnF 4 , we observe much more pronounced differences between MLP and DIPPIM. Firstly, the RMSE for the DIPPIM is nearly 14 times higher compared to that of the MLP. This higher RMSE for BaSnF 4 compared to NaF results from the presence of more complex local environments in BaSnF 4 , necessitating a more sophisticated and adaptable approach, which the MLP effectively addresses. Secondly, we observe that the distribution of errors for the DIPPIM presents some systematic errors, with a strong deviation from a Gaussian behavior, which means that some important interactions are missed during the fitting process. 11 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. However, for BaSnF 4 , we observe more significant differences as shown in FIG. 4. Although the position of the first peaks, and thus the interatomic distances, are well caught by the DIPPIM, this model is not able to reproduce well the shape in various cases. Notably, in the Sn-Sn RDF around 4 Å, the DIPPIM fails to capture the peak splitting. We attribute this discrepancy to the presence of stereoactive Sn lone pairs. The intensity and the shape 12 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

PLEASE CITE THIS ARTICLE AS

C. Vibrational properties

In the following section, we further compared the capability of MLPs and DIPPIMs to predict the vibrational properties of the materials by calculating the power spectrum and the phonon dispersion relations.

Power spectrum can be used to determine vibrational-rotational normal modes, or the vibrational density of states (VDOS) probed by experimental techniques such as infrared and Raman spectroscopies [START_REF] Wendler | Short time dynamics of ionic liquids in AIMD-based power spectra[END_REF] As for the structure, the differences become more significant in the case of BaSnF 4 . Unlike MLP potential, which is still capable of reproducing the peak positions and intensities, the DIPPIM almost fails to reproduce the DFT results. Specifically, in the Sn power spectra, there are two major peaks below and above 100 cm -1 . However, in the spectra obtained from DIPPIM, the peak above 100 cm -1 is almost missing and the peak below 100 cm -1 also shifts towards higher frequency. The failure shows the deficiency of DIPPIM to model the vibrational properties for BaSnF 4 .

14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. To go further, the phonon description is critical to understand the vibrational dynamics of a material and closely related to various thermal transport properties such as the thermal conductivity. We calculated the phonon dispersion of NaF with the same supercell as we did for our MLP training and DIPPIM fitting. The phonon dispersion results shown in FIG. 7 utilized the finite displacement method as implemented in Phonopy code [START_REF] Togo | First-principles Phonon Calculations with Phonopy and Phono3py[END_REF][START_REF] Togo | Implementation strategies in phonopy and phono3py[END_REF] . We can see that both methods agree well with the result from DFT.

PLEASE CITE THIS ARTICLE AS

On the other hand, the phonon dispersion undoubtedly demonstrates the superiority of MLP when applied to BaSnF 4 as shown in FIG. 8 . While the MLP phonon dispersion does not perfectly match that of DFT, it qualitatively captures the essential features. As for NaF, we calculated the phonon dispersion using a supercell of the same size to reproduce the DFT results. As we can see the phonon frequency from DIPPIM significantly differs from that of DFT. Note that the existence of imaginary modes in BaSnF 4 indicates that this structure is unstable at 0 K. This is consistent with experiments since it has been shown that the tetragonal phase of BaSnF 4 occurs upon annealing at 460 K 4 . 15 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

D. Equation of state

Finally, in order to test the transferability of the potentials within a range of physical conditions not encountered during the training/fitting procedure, we computed the 0 K equation of state (EOS) of NaF and BaSnF 4 . The EOS describes the material behavior under varied pressure and are of critical importance to the development and application of materials [START_REF] Latimer | Evaluation of thermodynamic equations of state across chemistry and structure in the materials project[END_REF] . As shown in FIG. 9, we strained the volume from -10% ∼ 10% and plotted the normalized energies with respect to the volume. Next, we fitted the energy-volume relationship to Birch-Murnaghan equation of state [START_REF] Birch | Finite Elastic Strain of Cubic Crystals[END_REF][START_REF] Murnaghan | The Compressibility of Media under Extreme Pressures[END_REF] and the fitted modulus of NaF is 44.70 GPa, 55.30 GPa and 55.83 GPa respectively for DFT, MLP and DIPPIM. In the case of BaSnF 4 , we obtained the following bulk modulus for DFT, MLP and DIPPIM: 68.89 GPa, 62.61 GPa and 115.58 GPa. The MLP computed bulk modulus deviates by 23.7% and 9.1% for NaF and BaSnF 4 from the DFT computed one, which can be ascribed to the fact that MLP has not been trained at various volumes. When pressure is close to zero, we have a nearly perfect agreement, which is reasonable because the thermodynamic condition is the same. For DIPPIM, the deviations of 24.8% and 67.8% are both larger than that of MLP, albeit the difference is almost negligible in the case of NaF. It is worth noting that the DIPPIM computed bulk modulus is astonishingly larger for BaSnF 4 , which again confirmed the failure of this potential in correctly describing potential energy surface of BaSnF 4 .

16

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. performed with MetalWalls, compiled on CPU as well. We can see that in general, the DIPPIM is more than 10 times more efficient than MLP. Therefore, from the computational point of view, DIPPIM remains the solution of choice for systems for which it is accurate, i.e. in the absence of specific interactions, while for more complex materials such as BaSnF 4 , the accuracy gain makes MLP a better choice than classical IP. 17 To further compare the two methods, we computed a series of structural (RDF), vibrational (power spectrum, phonon dispersion) and thermodynamic (EOS) observables. In all cases, the results demonstrate that our trained MLPs for both are robust and versatile. We also conclude that for simple systems like NaF, classical IP is sufficient to model the system and give accurate representation of the potential energy surface at high efficiency. While for materials with more sophisticated interactions as in BaSnF 4 , it is evident that fitting the interaction potential fails to provide a correct description of the system and using MLP is more advantageous.

PLEASE CITE THIS ARTICLE AS

SUPPLEMENTARY MATERIAL

The supplementary material contains the structural parameters obtained from geometry optimization using different functionals, the DIPPIM parameters for NaF and BaSnF 4 and a summary of the convergence test for the DIPPIM fitting for BaSnF 4 .

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

FIG. 1 :

 1 FIG. 1: (a) Crystal structure of BaSnF 4 (b) Crystal structure of NaF

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

( 2 )

 2 i (r), and number of radial basis functions for angular descriptors ρ

( 3 )

 3 i (r). The indices n, l, m represent the radial, angular and magnetic quantum numbers. The descriptors are expressed as coefficients vectors in the MLP:

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

( 2 )

 2 i (r) and angular descriptors ρ

( 3 )

 3 i (r) were set to 6 Å with a Gaussian broadening width of 0.5 Å. The number of spherical Bessel functions for radial (N 0 R ) and angular descriptors (N l R ) was set to 8 for both. L max , the maximum value of the angular quantum number, was set to 4. The kernel function sharpness parameter ζ[START_REF] Liu | Rechargeable anion-shuttle batteries for low-cost energy storage[END_REF] was set to 4. Ultimately, the algorithm sampled 244 and 516 reference structure datasets (from which the local reference configurations are extracted) for NaF and BaSnF 4 respectively. For NaF, 240 and 221 local reference configurations are selected for Na and F, indicating that Na and F local environments are equally complex. In the case of BaSnF 4 , 130, 606, and 1388 local reference configurations for Ba, Sn, and F have been sampled by the code. Intuitively, the Ba atoms have the most symmetric environment, thus requiring a smaller amount of reference configurations to describe them. In the case of Sn atoms, the presence of stereoactive lone pairs pointing towards the Sn-Sn empty layer results in the need for more reference configurations. Finally, it is not surprising that we need the highest number of local reference configurations for fluoride ions since they sit in various crystallographic site with different cationic environment. In addition they are the conducting species inside the structure and 7 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343 are highly mobile even at room temperature.

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

  DOI: 10.1063/5.0169343 Accepted to J. Chem. Phys. 10.1063/5.0169343

  FIG. 2: (a)NaF forces comparison between DFT and DIPPIM. (b)NaF forces comparison between DFT and MLP. (c)BaSnF 4 forces comparison between DFT and DIPPIM. (d) BaSnF 4 forces comparison between DFT and MLP.

  FIG. 3: Partial RDFs obtained at 298 K for NaF with AIMD (cyan, solid), DIPPIM (purple, dashed) and MLP (orange, dotted)

  FIG. 4: Partial RDFs obtained at 298 K for BaSnF 4 with AIMD (cyan, solid), DIPPIM (purple, dashed) and MLP (orange, dotted)

  FIG.5and FIG.6show the simulated mass-weighted power spectrum for NaF and BaSnF 4 , respectively. The power spectrum is calculated based on a separate 20 ps trajectory to ensure the convergence. For NaF, both DIPPIM and MLP can capture the normal nodes

  FIG. 7: Phonon dispersion of NaF obtained using DFT (cyan, solid), DIPPIM (purple, dashed) and MLP (orange, dotted)-based calculations

  FIG. 8: (a) Phonon dispersion of BaSnF 4 obtained using DFT (cyan, solid) and MLP (orange, dotted)-based calculations. (b)Phonon dispersion of BaSnF 4 obtained using DFT (cyan, solid) and DIPPIM (purple, dashed)-based calculations.

  FIG. 9: (a) EOS of NaF as obtained with DFT (cyan, circle), MLP (orange, triangle), and DIPPIM (purple, cross). (b) EOS of BaSnF 4 as obtained with DFT (cyan, circle),DIPPIM (purple, cross and MLP (orange, triangle))

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI: 10.1063/5.0169343IV. SUMMARYIn summary, we have trained MLP and fitted DIPPIM for both NaF and BaSnF 4 . In all the cases, the parameterization was made based on ab initio calculations only, i.e. no experimental information was used during this process. As a first accuracy test, we compared the forces computed from MLP and DIPPIM to those from DFT on the same set of configurations. Although the error made by the two models is similar in the case of NaF, the distribution of errors has a much more Gaussian behavior in the case of MLP, which tends to show that the DIPPIM fails to capture some of the interactions. The difference is much more pronounced in the case of BaSnF 4 due to the lack of a term to represent the stereoactive lone pairs of Sn in the DIPPIM.

TABLE I :

 I Computational cost for typical simulations performed with the MLP and DIPPIM It is worthwhile to mention the efficiency of the two interaction potentials. Typical timings are shown in Table. I. The calculations are performed on French supercomputer Jean-Zay where each node comprises of 40 cores 38 . The current efficiency tests for MLP were performed with the CPU compiled 6.3.2 version of VASP although the newer version is claimed to have significantly boosted efficiency. The efficiency tests for DIPPIM were

		number of atoms	number of nodes	MLP (s/step)	DIPPIM (s/step)
	NaF	1728	2	0.4907	0.0448
	BaSnF 4	600	3	0.3782	0.0285
	E. Computational performances

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI: 10.1063/5.0169343

ACKNOWLEDGMENTS

This work was granted access to the HPC resources of CINES under Allocation A0140910463 made by GENCI.

DATA AVAILABILITY

The data that support the findings of this study (input files for the simulations, raw data used for the various figures) are available on Zenodo (http://dx.doi.org/10.5281/zenodo.8354419).


published or not. The documents may come

A construction of the Shephard-Todd group G_32

through the Weyl group of type E_6

Cédric Bonnafé

A CONSTRUCTION OF THE SHEPHARD-TODD GROUP G 32

THROUGH THE WEYL GROUP OF TYPE E 6 by

CÉDRIC BONNAFÉ

Abstract. -It is well-known that the quotient of the derived subgroup of the Shephard-Todd complex reflection group G32 (which has rank 4) by its center is isomorphic to the derived subgroup of the Weyl group of type E6. We show that this isomorphism can be realized through the second exterior power, and take the opportunity to propose an alternative construction of the group G32.

Let G 32 denote the complex reflection group constructed by Shephard-Todd [ShTo] and let E 6 be a Weyl group of type E 6 (which is denoted by G 35 in the Shephard-Todd classification). Let Sp 4 (F 3 ) (resp. SO 5 (F 3 )) denote the symplectic (resp. orthogonal) group of dimension 4 (resp. 5) over the finite field with three éléments F 3 . Let Ω 5 (F 3 ) be the image of Sp 4 (F 3 ) in SO 5 (F 3 ) through the natural morphism Sp 4 (F 3 ) → SO 5 (F 3 ): this is the normal subgroup of index 2 of SO 5 (F 3 ). Finally, if G is a group, let D(G) and Z(G) denote respectively its derived subgroup and its center and, if d is a non-zero natural number, let µ d be the group of d-th roots of unity in C.

It is shown in [START_REF] Lehrer | Unitary reflection groups[END_REF]théo. 8.43 et 8.54] that G 32 ≃ µ 3 × Sp 4 (F 3 ) and that E 6 ≃ SO 5 (F 3 ). In particular,

( * ) D(G 32 )/µ 2 ≃ Ω 5 (F 3 ) ≃ D(E 6 ).
The purpose of this note is to present a direct elementary explanation of the isomorphism D(G 32 )/µ 2 ≃ D(E 6 ), which in fact allows us to construct the the complex reflection group G 32 from the rational reflection group E 6 . This construction uses the classical morphism SL 4 (C) -→ SO 6 (C), and follows the same lines as our previous paper [Bon] (in which we constructed the complex reflection group G 31 from the Weyl group of type B 6 ). This note does not pretend to prove a deep result: it is just a a nice example of the application of the classical theory of reflection groups (invariants, Springer theory,...).

Remark.

It is shown in [START_REF] Lehrer | Unitary reflection groups[END_REF]Theo. 8.53] that G 33 ≃ µ 2 × Ω 5 (F 3 ). In particular, this gives an indirect isomorphism D(G 33 ) ≃ D(E 6 ). However, we do not know of any construction of this isomorphism which would be in the same spirit as above.

The morphism SL

4 (C) → SO 6 (C)
We recall here the construction of this morphism, using some notation from [Bon]. Let us fix a complex vector space V of dimension 4 and let

Λ : GL(V ) -→ GL(∧ 2 V ) g -→ ∧ 2 g
be the natural morphism of algebraic groups. Note that ∧ 2 V has dimension 6 and that

(1.1)

Ker Λ ≃ µ 2 = {± Id V }.
The lists of eigenvalues of elements of GL(V ) or GL(∧ 2 V ) will always be given with multiplicities. If g ∈ GL(V ) admits a, b, c and d as eigenvalues, then

(1.2) Λ(g) admits ab, ac, ad, bc, bd and cd as eigenvalues.

In particular,

(1.3) det Λ(g) = (det g) 3 .
Let us fix now a generator ε of the one-dimensional vector space ∧ 4 V . The choice of this generator allows to identify C and ∧ 4 V and to define a bilinear form

β ∧ : ∧ 2 V × ∧ 2 V -→ C (x, y) -→ x ∧ y.
This bilinear form is symmetric and non-degenerate. By definition of the determinant,

(1.4) β ∧ (Λ(g)(x), Λ(g)(y)) = (det g)β ∧ (x, y)
for all g ∈ GL(V ) and x, y ∈ ∧ 2 V . For dimension and connectedness reasons, the image of

GL(V ) through Λ is the neutral component CO(∧ 2 V ) • of the conformal orthogonal group CO(∧ 2 V ) = CO(∧ 2 V, β ∧ )
and Λ induces an isomorphism of algebraic groups

(1.5) SL(V )/µ 2 ≃ SO(∧ 2 V ).
The next elementary lemma will be useful (in this paper, we denote by j a primitive third root of unity).

Lemma 1.6. -Let g ∈ SO(∧ 2 V ) having j, j, j, j 2 , j 2 and j 2 as eigenvalues. Then there exists a unique element g ∈ SL(V ), having 1, j, j and j as eigenvalues, and such that

Λ -1 ({g, g -1 }) = {±g, ±g -1 }. Proof. -Let h ∈ SL(V ) be such that Λ(h) = g. Then Λ -1 ({g, g -1 }) = {±h, ±h -1 }. Let a,
b, c and d be the eigenvalues of h. By (1.2), by reordering if necessary the eigenvalues of h, we may assume that ab = ac = j. In particular, b = c. So bd = cd, which implies that bd = cd = j 2 (for otherwise j would be an eigenvalue of g with multiplicity ⩾ 4). Still by (1.2), (ad, bc) = (j, j 2 ) or (ad, bc) = (j 2 , j).

If (ad, bc) = (j, j 2 ), then b = c = d and b 2 = j 2 , so b = c = d = ηj and a = η, for some η ∈ {±1}. By replacing h by -h if necessary, the eigenvalues of h are then 1, j, j, j and so h is indeed the unique element of Λ -1 ({g, g -1 }) admitting this list of eigenvalues.

If (ad, bc) = (j 2 , j), then d = jb = jc and b 2 = j, which implies that a = b = c = ηj 2 and d = η, for some η ∈ {±1}. By replacing h by -h if necessary, the eigenvalues of h are then 1, j 2 , j 2 , j 2 and so h -1 is indeed the unique element of Λ -1 ({g, g -1 }) admitting 1, j, j, j as list of eigenvalues.

Corollary 1.7. -Let g ∈ SO(∧ 2 V ) having j, j, j, j 2 , j 2 and j 2 as eigenvalues. Then g and g -1 are not conjugate in SO(∧ 2 V ).

Proof. -Assume that we have found w ∈ SO(∧ 2 V ) such that wgw -1 = g -1 . Let ω and γ be respective preimages of w and g in SL(V ). Then Λ(ωγω -1 ) = g -1 = Λ(γ -1 ). This shows that γ -1 is conjugate to γ or to -γ, which is impossible by examining the possible lists of eigenvalues of γ obtained in the proof of Lemma 1.6.

Construction of G 32

Let us see E 6 as a finite subgroup of O(∧ 2 V ) generated by reflections. Set

jSO(∧ 2 V ) = ⟨SO(∧ 2 V ), j Id ∧ 2 V ⟩ and 3 √ SL(V ) = {g ∈ GL(V ) | det(g) 3 = 1} = ⟨SL(V ), j Id V ⟩. Then (2.1) Λ -1 jSO(∧ 2 V ) = 3 √ SL(V ).
We then define

W = Λ -1 (⟨D(E 6 ), j Id ∧ 2 V ⟩).
It is a subgroup of 3 √ SL(V ). The aim of this note is to show that W is isomorphic to the complex reflection group G 32 of Shephard-Todd. Note first that (2.2) µ 6 ⊂ W.

2.A. Reflections in W .

-The list of degrees E 6 is 2, 5, 6, 8, 9, 12 while its list of codegrees is 0, 3, 4, 6, 7, 10 (see [START_REF] Broué | Introduction to complex reflection groups and their braid groups[END_REF]Table A.3]). In particular, exactly 3 of the degrees are divisible by 3, which shows [START_REF] Springer | Regular elements of finite reflection groups[END_REF]Theo. 3.4] that E 6 contains an element w 3 admitting the eigenvalue j with multiplicity 3. We will denote by C 3 the conjugacy class of w 3 in W . Since also exactly 3 of the codegrees are divisible by 3, this implies, for instance by [START_REF] Lehrer | Invariant theory and eigenspaces for unitary reflection groups[END_REF]Theo. 1.2], that w 3 is regular in the sense of Springer [START_REF] Springer | Regular elements of finite reflection groups[END_REF]§4] (that is, admits an eigenvector for the eigenvalue j whose stabilizer in E 6 is trivial). Since E 6 is a rational group, w 3 also admits j 2 as an eigenvalue with multiplicity 3. Hence, the eigenvalues of w 3 are j, j, j, j 2 , j 2 and j 2 . In particular, det w 3 = 1 and so

w 3 ∈ D(E 6 ) = E 6 ∩ SO(∧ 2 V ).
By [START_REF] Springer | Regular elements of finite reflection groups[END_REF]Theo. 4.2(iii)], the centralizer of w 3 in E 6 has order 6 • 9 • 12, which shows that

(2.3) |C 3 | = 2 • 5 • 8 = 80.
Moreover [START_REF] Springer | Regular elements of finite reflection groups[END_REF]Theo. 4.2(iii)], if w ∈ E 6 , then

(2.4) w ∈ C 3 if and only if dim Ker(w -j Id V ) = 3.

Hence, if w ∈ C 3 , then w -1 ∈ C 3 but w -1 is not conjugate to w in D(E 6 ) (by Corollary 1.7). Now, let Ref(W ) denote the set of reflections of W . If s ∈ Ref(W ), the eigenvalues Λ(s) are 1, 1, 1, det(s), det(s), det(s). Since det(s) ∈ {j, j 2 }, the eigenvalues of det(s)Λ(s) are then j, j, j, j 2 , j 2 , j 2 . It follows from (2.4) that det(s)Λ(s) ∈ C 3 . This defines a map

λ : Ref(W ) -→ C 3 s -→ det(s)Λ(s).
The next result will be very useful:

Lemma 2.5. -The map λ is bijective.

Proof. -First, if λ(s) = λ(s ′ ), then there exists ξ ∈ C × such that s ′ = ξs. Since s and s ′ are reflections, this is possible only if ξ = 1, and so s = s ′ . This shows that λ is injective.

Let us now show the surjectivity. Let w ∈ C 3 . By Corollary 1.7, there exists a unique w ∈ Λ -1 ({w, w -1 }) admitting 1, j, j, j as eigenvalues. Then j 2 w and j w-1 are reflections satisfying λ(j 2 w) = det(j 2 w)Λ(j 2 w) = j 8 • j 4 Λ( w) = Λ( w) and λ(j w-1 ) = det(j w-1 )Λ(j w-1 ) = j 4 • j 2 Λ( w-1 ) = Λ( w) -1 .

So w ∈ {λ(j 2 w), λ(j w-1 )}, which shows that λ is surjective.

We then deduce from (2.3) and Lemma 2.5 that

(2.6) |Ref(W )| = 80.

2.B.

Structure of W . -Our main result is Theorem 2.7 below: the proof we propose here uses neither known properties of the group G 32 nor the classification of complex reflection groups and so might be viewed as an alternative construction of G 32 starting from E 6 (however, note that we use properties of E 6 ).

Theorem 2.7. -The group W : (a) has order 155 520; (b) is generated by reflections of order 3; (c) is irreducible and primitive;

(d) admits 12, 18, 24, 30 as list of degrees.

Proof. (2.9)

-Set E # 6 = ⟨D(E 6 ), j Id ∧ 2 V ⟩ and W + = W ∩ SL(V ). Recall that |E 6 | = 51 840. So (2.8) |D(E 6 )| =
Z(W ) = µ 6 and W/µ 6 ≃ D(E 6 ). Let R = {det(s) -1 s | s ∈ Ref(W )}. Set G = ⟨Ref(W )⟩ and H = ⟨R⟩.
The statement (b) is equivalent to the following one

(#) W = G.
First, Λ(R) = C 3 and so Λ(H) = D(E 6 ) (as this last group is simple and C 3 is a conjugacy class). In particular, W = H • µ 6 and so, since H ⊂ G • µ 3 , we obtain

W = G • µ 6 ,
which is almost the expected result (#). Before showing (#), note that, since Λ(H) = D(E 6 ), we have that H (and so G) acts irreducibly on V . Moreover, if G is not primitive, then G (and so H) would be monomial [START_REF] Lehrer | Unitary reflection groups[END_REF]Lem. 2.12], which would imply that Λ(H) = D(E 6 ) is monomial, which is false. So From the second equality, we deduce that at least one of the e i 's (say e 3 ) is odd. From the first equality, we deduce that at least one of the e i 's (say e 2 ) is even and so e 1 is also even.

The first equality shows that e 1 or e 2 (say e 2 ) is divisible by 8. So e 2 ∈ {8, 16, 24}. A quick inspection of the possibilities shows that e 1 = 4, e 2 = 16 and {e 3 , e 4 } = {3, 5}. The degrees of G are then 9, 12, 15 and 48. Since 16 divides one of the degrees, it then follows from [START_REF] Springer | Regular elements of finite reflection groups[END_REF]théo. 3.4(i)] that G = H × µ 3 contains an element of order 16 and so H contains an element of order 16. Therefore, D(E 6 ) = Λ(H) contains an element of order 8, which is impossible (see Remark 2.11 below for a proof of this fact based only on Springer theory). This contradicts the fact that |G| = 77 760. So we have shown (#), that is,

W = ⟨Ref(W )⟩ = G.
This is statement (b). In particular, µ 6 ⊂ W and so all the d i 's are divisibles by 6. Set

a i = d i /6. Then (♢) implies that a 1 a 2 a 3 a 4 = 120 = 2 3 • 3 • 5, a 1 + a 2 + a 3 + a 4 = 14,
By the same argument as before, since µ 12 ̸ ⊂ W , we may assume that a 3 and a 4 are odd and that a 1 and a 2 are even. A quick inspection of the possibilities shows that {a 1 , a 2 } = {2, 4} and {e 3 , e 4 } = {3, 5}. This concludes the proof of (d).

Corollary 2.10. -The group W is isomorphic to the reflection group G 32 of Shephard-Todd.

Proof. -This follows from Theorem 2.7 and from the classification of complex reflection groups [ShTo].

Hence, Corollary 2.10 gives an explanation for the fact, mentioned in the introduction, that D(G 32 )/µ 2 ≃ D(E 6 ): the isomorphism is realized by Λ.

Remarque 2.11. -In the proof of Theorem 2.7, we have used the fact that D(E 6 ) does not contain any element of order 8. This fact can be easily obtained by a computer calculation for instance, but we propose here a proof using only Springer theory. Let w ∈ E 6 be an element of order 8. Then w necessarily admits an eigenvalue which is a primitive 8-th root of unity ζ. As only one of the degrees of E 6 and only one of the codegrees of E 6 is divisible by 8, this implies that w is a regular element in the sense of Springer [START_REF] Lehrer | Invariant theory and eigenspaces for unitary reflection groups[END_REF]Theo. 1.2]. Then [START_REF] Springer | Regular elements of finite reflection groups[END_REF]Theo. 4.2(v)] the list of eigenvalues of w is ζ -1 , ζ -4 , ζ -5 , ζ -7 , ζ -8 , ζ -11 , and so det(w) = ζ -36 = -1. So w ∈ D(E 6 ).

  25 920, |E # 6 | = 77 760, |W | = 155 520 and |W + | = 51 840. This shows (a). Moreover,

  this proof by showing simultaneously (b) et(d). Let d 1 , d 2 , d 3 , d 4 be the degrees of G. It follows from (2.6) and for instance from[START_REF] Broué | Introduction to complex reflection groups and their braid groups[END_REF] Theo. 4.1] that(♢) d 1 d 2 d 3 d 4 = |G|, d 1 + d 2 + d 3 + d 4 = 84.The morphism det : G → µ 3 is surjective (since det(s) ∈ {j, j 2 } for any s ∈ Ref(W )), and this implies thatµ 3 ⊂ G (because G/(G ∩ µ 6 ) ≃ D(E 6 ) is simple). Itremains to show that |G| ̸ = |W |/2 = 77 760. So assume that |G| = 77 760. Since µ 3 ⊂ G, all the d i 's are divisible by 3 and, since µ 6 ̸ ⊂ G, at least one of them (say d 4 ) is not divisible by 6. Write e i = d i /3. 2 e 3 e 4 = 960 = 2 6 • 3 • 5, e 1 + e 2 + e 3 + e 4 = 28, e 4 is odd.


D V Louzguine-Luzgin

K Georgarakis

J Andrieux

L Hennet

T Morishita

K Nishio

R V Belosludov

An atomistic study of the structural changes in a Zr-Cu-Ni-Al glass-forming liquid on vitrification monitored in-situ by X-ray diffraction and molecular dynamics simulation

Keywords: metallic glass, liquid, fragility, X-ray diffraction, structure, synchrotron radiation

come

An atomistic study of the structural changes in a Zr-Cu-Ni-Al glass-forming liquid on vitrification monitored in-situ by X-ray diffraction and molecular dynamics simulation

Introduction

The high glass-forming ability (GFA) of some alloy compositions allowed production of bulk metallic glasses (BMGs) (also called glassy alloys) in the thickness range of 1-100 mm using various casting processes [START_REF] Suryanarayana | Bulk metallic glasses[END_REF][START_REF] Louzguine-Luzgin | Metallic glasses and those composites[END_REF]. Zr-based BMGs are among the best metallic glass-formers known to date [START_REF] Suryanarayana | Bulk metallic glasses[END_REF][START_REF] Louzguine-Luzgin | Metallic glasses and those composites[END_REF][START_REF] Petrzhik | On conditions of bulk and surface glass formation of metallic alloys[END_REF] which thermal stability is also good [START_REF] Louzguine-Luzgin | On long-term stability of metallic glasses[END_REF]. Formation of a glassy phase in both Zr-Cu and Zr-Ni alloys is strongly enhanced by the addition of Al thus forming ternary and quaternary bulk metallic glasses with good glass-forming ability [START_REF] Suryanarayana | Bulk metallic glasses[END_REF][START_REF] Louzguine-Luzgin | Metallic glasses and those composites[END_REF][START_REF] Petrzhik | On conditions of bulk and surface glass formation of metallic alloys[END_REF][START_REF] Louzguine-Luzgin | On long-term stability of metallic glasses[END_REF][START_REF] Molokanov | Glass forming ability, structure and properties of Ti and Zrintermetallic compound based alloys[END_REF]. The atomic structure of these BMGs at room temperature was studied in detail [START_REF] Sato | Effect of Al on local structures of Zr-Ni and Zr-Cu metallic glasses[END_REF][START_REF] Antonowicz | Atomic structure of Zr-Cu-Al and Zr-Ni-Al amorphous alloys[END_REF].

Liquids retain their volume as crystals but flow under the action of gravity [ 8 , 9 ]. They constantly undergo a restructuring of the atomic structure [START_REF] Angell | Relaxation in liquids, polymers and plastic crystals-strong/fragile patterns and problems[END_REF][START_REF] Egami | Real-space description of dynamics of liquids[END_REF]. Above the equilibrium liquidus temperature (T l ) and slightly below it metallic melts show nearly Arrhenius-type temperature dependence of viscosity ():

) / exp( 0 RT E a    (1)
where  0 is a pre-exponential factor, R is the gas constant and E a is an activation energy [START_REF] Gale | Smithells Metals Reference Book 8-th[END_REF]. However, on cooling below the liquidus temperature, starting from a crossover temperature [START_REF] Blodgett | Proposal for universality in the viscosity of metallic liquids[END_REF], they exhibit a non-Arrhenius temperature dependence on viscosity which is known as fragility of liquids [START_REF] Angell | Relaxation in liquids, polymers and plastic crystals -strong/fragile patterns and problems[END_REF][START_REF] Hu | Liquid fragility. A key to going deep into materials of glassy states[END_REF].

Here one should mention that anomalous variation in a liquid viscosity in the vicinity of T l and above it was observed for Fe- [START_REF] Popel | Microheterogeneity of liquid metallic solutions and its influence on the structure and properties of rapidly quenched alloys[END_REF][START_REF] Dahlborg | Structure and properties of some glass-forming liquid alloys[END_REF] and Zr-based glass-forming alloys [START_REF] Polukhin | Presolidification changes in the structural-dynamic characteristics of glass-forming metallic melts during deep cooling, vitrification, and hydrogenation[END_REF] as well as for Al-based crystalline alloys [START_REF] Popel | Conditions for the formation of metastable phases during crystallization of Al-Zr alloys[END_REF]. The authors stated that liquid metallic alloys can have various structural states and significant overheating is required to dissolve clusters inherited from the solid state and to reach equilibrium liquid structure. The existence of liquid-liquid transitions in deeply supercooled state found for water [START_REF] Mishima | The relationship between liquid, supercooled and glassy water[END_REF], molecular liquids [START_REF] Katayama | A firstorder liquid-liquid phase transition in phosphorus[END_REF][START_REF] Kobayashi | The reversibility and first-order nature of liquid-liquid transition in a molecular liquid[END_REF] (confirmed by computer simulation [START_REF] Morishita | Liquid-liquid phase transitions of phosphorus via constant-pressure first-principles molecular dynamics simulations[END_REF]) and even metallic glasses [START_REF] Stolpe | Structural changes during a liquid-liquid transition in the deeply undercooled Zr 58.5 Cu 15.6 Ni 12.8 Al 10.3 Nb 2.8 bulk metallic glass forming melt[END_REF] should also be taken into account when describing fragility.

The inverse temperature -logarithmic viscosity plot [START_REF] Laughlin | Viscous flow in simple organic liquids[END_REF][START_REF] Angell | Formation of glasses from liquids and biopolymers[END_REF] with some limitations [START_REF] Louzguine-Luzgin | On limitations of the viscosity versus temperature plot for glass-forming substances[END_REF][START_REF] Ikeda | Analysis and characterization of the transition from the Arrhenius to non-Arrhenius structural relaxation in fragile glass-forming liquids[END_REF] is a very useful illustration of the difference between so-called "strong" and "fragile" liquids [ 29 , 30 ].

Fragility is also related to the glass-forming ability of BMGs [START_REF] Senkov | Correlation between fragility and glass-forming ability of metallic alloys[END_REF][START_REF] Louzguine-Luzgin | Role of different factors in the glass-forming ability of binary alloys[END_REF] (with some argumentation [START_REF] Johnson | Quantifying the origin of metallic glass formation[END_REF]) partly because stronger liquids, in general, have a higher viscosity in the entire temperature range from T l down to T g . There are also experimental results which indicate correlation between liquid fragility and vibrational properties of the glass below T g [START_REF] Scopigno | Is the fragility of a liquid embedded in the properties of its glass?[END_REF]. The fragility index (m) [START_REF] Johnson | Quantifying the origin of metallic glass formation[END_REF][START_REF] Kim | Manipulation of thermal and mechanical stability by addition of multiple equiatomic rare-earth elements in Al-TM-RE metallic glasses[END_REF] of the supercooled liquid is calculated slightly above the glass-transition temperature (T g ) as a derivative:

m=dlog()/d(T g /T) (2).
Another way of representation of fragility is a famous Vogel-Fulcher-Tammann-Hesse equation [START_REF] Suryanarayana | Bulk metallic glasses[END_REF]:

= 0 exp[D*T 0 /(T-T 0 )] ( 3
)
where  0 , D (as an indicator of the fragility) and T 0 are the fitting parameters. This equation is widely used but describes the temperature dependence of viscosity well only in an intermediate temperature region [START_REF] Kawamura | Newtonian viscosity of supercooled liquid in a Pd 40 Ni 40 P 20 metallic glass[END_REF]. While some other equations make a better representation of the entire viscosity plot [START_REF] Cohen | Liquid-glass transition, a free-volume approach[END_REF][START_REF] Ojovan | Ordering and structural changes at the glass-liquid transition[END_REF] they have a larger number of fitting parameters which physical meaning could be less clear.

On the other hand fragility can be expressed by the ratio of the activation energies for viscous flow in the equation (1) at low (E L ), (close to T g ) and high (E H ) temperature region (slightly above T l ) [START_REF] Ojovan | Viscous flow and the viscosity of melts and glasses[END_REF][START_REF] Sanditov | Relaxation aspects of the liquid-glass transition[END_REF]:

R D = E L /E H (4).
Other methods are also used to describe and separate thermodynamic and kinetic fragility [START_REF] Kelton | Kinetic and structural fragility -A correlation between structures and dynamics in metallic liquids and glasses[END_REF].

Metallic glasses exhibit structural changes in the supercooled liquid state and such changes are found to be responsible for the liquid fragility. A relatively fragile Pd 42.5 Cu 30 Ni 7.5 P 20 melt (its fragility index m is close to 60 [START_REF] Fan | Thermodynamics, enthalpy relaxation and fragility of the bulk metallic glass-forming liquid Pd 43 Ni 10 Cu 27 P 20[END_REF]) during cooling was studied by using the real-space pair distribution function (G) (which is a reduced radial distribution function (RDF)) [START_REF] Louzguine-Luzgin | Comparative analysis of the structure of palladium-based bulk metallic glasses prepared by treatment of melts with flux[END_REF] and compared to the results obtained by the first-principles calculations. As a result strong correlation between the structural changes and the variation of viscosity was observed in the supercooled liquid [START_REF] Louzguine-Luzgin | Structural basis for supercooled liquid fragility established by synchrotronradiation method and computer simulation[END_REF]. The rate of structural rearrangements on cooling was enhanced in the supercooled liquid close to T g . It was also demonstrated for Ni-Nb-Ta [START_REF] Mauro | Structural evolution in Ni-Nb and Ni-Nb-Ta liquids and glasses -A measure of liquid fragility?[END_REF] and other [START_REF] Mauro | A structural signature of liquid fragility[END_REF][START_REF] Gangopadhyay | Recent progress in understanding high temperature dynamical properties and fragility in metallic liquids, and their connection with atomic structure[END_REF] liquid alloys.

Structural changes, thermal expansion and volume changes (structural relaxation) in Zr-based metallic glasses upon heating [START_REF] Yavari | Excess free volume in metallic glasses measured by X-ray diffraction[END_REF][START_REF] Mizuno | High-energy x-ray diffraction study of liquid structure of metallic glass-forming Zr 70 Cu 30 alloy[END_REF] and cooling [START_REF] Qu | Undercooling behavior of Zr-Cu-Ni-Al bulk metallic glasses investigated by in situ synchrotron high energy X-ray diffraction[END_REF] were studied by synchrotron X-ray radiation diffraction and transmission electron microscopy [START_REF] Louzguine | Formation of a nanoquasicrystalline phase in Zr-Cu-Ti-Ni metallic glass[END_REF][START_REF] Aronin | Nanocrystallization of bulk Zr-Cu-Ti metallic glass[END_REF]. Recently it was shown that even scanning tunneling microscopy allows achieving atomic scale resolution for metallic glasses [START_REF] Belosludov | The atomic structure of a bulk metallic glass resolved by scanning tunneling microscopy and ab-initio molecular dynamics simulation[END_REF].

Zr-Cu(Ni)-Al alloys are among the most studied bulk metallic glasses [START_REF] Louzguine-Luzgin | Exceptionally high nanoscale wear resistance of a Cu 47 Zr 45 Al 8 metallic glass with native and artificially grown oxide[END_REF]. A detailed study on vitrification of a ternary Zr 60 Cu 30 Al 10 alloy was performed recently in-situ by high energy synchrotron radiation X-ray diffraction from above the liquidus T l to room temperature [START_REF] Georgarakis | Probing the structure of a liquid metal during vitrification[END_REF]. Short (SRO) and medium range order (MRO) develop significantly during cooling the liquid phase to the glassy state.

Significant glassy structure and volume changes are found during severe plastic deformation [START_REF] Gunderov | High pressure torsion induced structural transformations in Ti-and Zr-based amorphous alloys[END_REF] and rolling [START_REF] Abrosimova | Influence of deformation on the structural transformation of the Pd 40 Ni 40 P 20 amorphous phase[END_REF]. Structural changes in Cu-Zr [START_REF] Mendelev | Experimental and computer simulation determination of the structural changes occurring through the liquid-glass transition in Cu-Zr alloys[END_REF][START_REF] Lekka | Clustering, microalloying and mechanical properties in Cu/Zr-based glassy models by molecular dynamics simulations and ab-initio computations[END_REF] and Cu-Zr-Al glasses [START_REF] Cheng | Alloying strongly influences the structure, dynamics, and glass forming ability of metallic supercooled liquids[END_REF] were also studied by MD computer simulation.

Here we study the atomic structure changes in a relatively strong Zr 55 Cu 30 Ni 5 Al 10 glass-forming liquid (one of the best glass-formers [START_REF] Yokoyama | Production of Zr 55 Cu 30 Ni 5 Al 10 glassy alloy rod of 30mm in diameter by a cap-cast technique[END_REF]) on vitrification and compare the results with those obtained for fragile metallic glasses. Its measured fragility parameter m for this liquid varies from about 45 [START_REF] Yamasaki | Viscosity measurements of Zr 55 Cu 30 Al 10 Ni 5 and Pd 40 Cu 30 Ni 10 P 20 supercooled liquid alloys by using a penetration viscometer[END_REF][START_REF] Qiao | Mechanical relaxation in a Zr-based bulk metallic glass: Analysis based on physical models[END_REF] to 29 [START_REF] Haruyama | Comparison of structural relaxation behavior in as-cast and pre-annealed Zr-based bulk metallic glasses just below glass transition[END_REF] but it is lower than that obtained for the Pd-Cu-Ni-P one studied earlier.

Experimental Procedure

An ingot of the Zr 55 Cu 30 Ni 5 Al 10 alloy (the composition is given in nominal atomic percentage) was prepared from pure metals of more than 99.9 mass% purity using the arc-melting method.

Thermal stability of the samples was tested using a differential scanning calorimeter (DSC).

In-situ X-ray diffraction experiments were carried out at the European Synchrotron Radiation Facility (ESRF). The incident beam wavelength was 0.1245 Å (100 keV). Diffraction spectra were acquired in transmission mode by a 2D detector. An aerodynamic levitation setup [START_REF] Hennet | Aerodynamic levitation and laser heating: applications at synchrotron and neutron sources[END_REF] was used to heat up, melt and cool down the metallic glass spheres in situ under synchrotron radiation. Prior to experiments, the levitation chamber was first evacuated down to 10 -5 mbar and then filled with the levitation gas. The spheres were levitating by high purity argon (containing less than 0.01 mass ppm O2) flowing through an aluminum nozzle and melted using two CO 2 laser beams. Two pyrometers were used to record the temperature while temperature values below 620 K were calculated.

Diffraction intensities were recorded during cooling with a time resolution of 100 ms. The radial intensities were integrated over 90 degrees (1/4 of the entire halo) to obtain one-dimensional spectrum and enhance the signal to noise ratio. The accuracy of the temperature values is estimated at ±10 K for the data above 900 K and better than ±7 K below 900 K.

After necessary corrections the measured intensity was converted to electron units per atom with the generalized Krogh-Moe-Norman method [START_REF] Waseda | The structure of non-crystalline materials[END_REF]. The density value of 6800 kg/m 3 for the alloy was obtained from Ref. [START_REF] Haruyama | Volume and enthalpy relaxation in Zr 55 Cu 30 Ni 5 Al 10 bulk metallic glass[END_REF]. The total structure factor S(Q) and the interference function Qi(Q) (Q = 4sinθ/λ, θ is the diffraction angle) were obtained from the coherent scattering intensity by using atomic scattering factors. The radial distribution RDF(R) and pair distribution functions PDF(R) and G(R) were obtained by the Fourier transform of Qi(Q) [START_REF] Waseda | The structure of non-crystalline materials[END_REF]. The maximum scattering vector (Q) value obtained was 250 nm -1 . However, as the noise level at high Q values was rather high the Fourier

Transform was performed only until 150 nm -1 . Termination of the Fourier Transform at higher Q max values was found to increase the noise of G(r) without any positive effect on the curves.

The atomic structure of the Zr 55 Cu 30 Ni 5 Al 10 alloy in both liquid and glassy states was also modeled using the first-principles molecular dynamics simulation. The calculations were performed using density functional theory based pseudopotential plane-wave method as implemented in VASP code [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. The atomic configuration was modeled within a cubic supercell consisting of 500 atoms.

The spin polarized generalized gradient corrected Perdew-Burke-Ernzerhof exchange-correlation functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and the projector augmented wave method [START_REF] Blöchl | Projector augmented-wave method[END_REF] were used in order to accurately describe the interaction between the ion and electron. The Brillouin zone integrations were carried out using the -point. At the first step, the random atomic structure was melted at 2000 K for a total of 10 ps. After that, the obtained system in liquid state was additionally equilibrated for 10 ps at 1400 K and 1000 K. The additional volume relaxations were performed for these structures using stress tensor calculation. In order to obtain the structures at lower temperatures owing to slower relaxation the cell obtained at 1000 K was cooled down to 800 K and then to 500 K with a cooling rate of 2×10 14 K/s.

The volume optimization was done for structures at 800 K and 500 K. After cell optimizations, all structures were equilibrated again for 2.5 ps at the same temperatures, respectively. The calculation of electron density of states were performed using enlarged Brillouin zone of 6 × 6 × 6 in order to get denser electron coverage.

Results and discussion

The alloy sample was heated up, melted and then cooled down at about 85 K/s (average cooling rate) from 1335 K (T l ~ 1150 K) to 388 K and the alloy vitrified [START_REF] Haruyama | Volume and enthalpy relaxation in Zr 55 Cu 30 Ni 5 Al 10 bulk metallic glass[END_REF]. The measured intensity and the structure factor profiles obtained at two temperatures are shown in Fig. 1a. The glasstransition temperature measured by DSC at the heating rate of 0.67 K/s was about 680 K. This value is close to those found in the earlier works [START_REF] Haruyama | Volume and enthalpy relaxation in Zr 55 Cu 30 Ni 5 Al 10 bulk metallic glass[END_REF][START_REF] Sergiienko | Investigation of Zr 55 Cu 30 Al 10 Ni 5 bulk amorphous alloy crystallization[END_REF].

PDF(R)=G(R)/(4 0 R)+1, ( 5
)
where  0 is the atomic number density, are also shown in Fig. 1c. In order to separate influence of different atomic pairs the G(R) peaks were fitted with the several Gaussian functions. Fitting of four G(R) maxima (the baseline was corrected to make both outermost points of each coordination shell maximum be equal to zero) using 7 Gaussian peak functions produced a reasonable correspondence with the original G(R) plot (Fig. 2a). The temperature dependence of the center of mass positions of 5 th and 6 th peaks (3 th coordination shell) indicates linear thermal contraction on cooling with a linear thermal expansion coefficient (a l ) of 2.810 -5 K -1 in the liquid and of 1.110 -5 K -1 in the glassy state (Fig. 2b). These values are quite close to those obtained by density measurements [START_REF] Haruyama | Volume and enthalpy relaxation in Zr 55 Cu 30 Ni 5 Al 10 bulk metallic glass[END_REF]. Opposite to the peaks related to other coordination shells (3 rd , 4 th , etc...) P1 and P2 of the first coordination shell shift to higher values on undercooling (Fig. 3a,b) as was found earlier for the Pd-Cu-Ni-P [START_REF] Mattern | Structural behavior of Pd 40 Cu 30 Ni 10 P 20 bulk metallic glass below and above the glass transition[END_REF][START_REF] Georgarakis | Variations in atomic structural features of a supercooled Pd-Ni-Cu-P glass forming liquid during in situ vitrification[END_REF], Pd-Si [START_REF] Louzguine-Luzgin | Atomic structure changes and phase transformation behavior in Pd-Si bulk glass-forming alloy[END_REF], Cu-Zr [START_REF] Mattern | Temperature dependence of the short-range order of Cu 65 Zr 35 metallic glass[END_REF], Zr-Cu-Ni-Al [START_REF] Tong | Structural evolution and strength change of a metallic glass at different temperatures[END_REF], Fe-B [START_REF] Louzguine-Luzgin | Structural changes in liquid Fe and Fe-B alloy on cooling[END_REF], many other alloys [START_REF] Gangopadhyay | Anomalous thermal contraction of the first coordination shell in metallic alloy liquids[END_REF] as well as for pure metals [START_REF] Lou | Negative expansions of interatomic distances in metallic melts[END_REF]. On cooling the present Zr-Cu-Ni-Al liquid alloy exhibits contraction in 3 rd (P5 and P6) and 4 th (P7) coordination shells (Fig. 2b). The peak position in the 2 nd coordination shell, in general, did not change with temperature. Although P3 and P4 move in opposite directions (Fig. 3c) the average center of mass value of P3 and P4 is about 0.54 nm in the entire temperature range. As it has been already mentioned before, seven Gaussian function peaks represent very well the entire shape of G(R) (Fig. 2a). However, in the earlier works it was shown that owing to a complex asymmetrical shape of the first peak the maximum of G(R) does not necessarily characterize all interatomic distances and the entire shape of the first peak should be taken into consideration [START_REF] Ding | Temperature effects on atomic pair distribution functions of melts[END_REF][START_REF] Sukhomlinov | Determination of accurate mean bond lengths from radial distribution functions[END_REF].

As one can see in Fig. 1b On subsequent heating of the glassy sample vitrified on cooling, the change of slope of temperature dependence of the first and second peaks (P1 and P2) takes place at about 500 K (see supplementary data, please) indicating densification of the glassy phase. This temperature corresponds well to the temperature of beginning of structural relaxation of the glassy phase on heating.

As in case of the Pd-Cu-Ni-P alloy studied earlier [START_REF] Louzguine-Luzgin | Structural basis for supercooled liquid fragility established by synchrotronradiation method and computer simulation[END_REF] AP1/AP2 (the area under P1 divided by the area under P2) ratio is nearly constant in the equilibrium liquid state. However, it changes on cooling in the supercooled liquid state (Fig. 3a). It indicates that the number of Zr-Cu and Zr-Ni atomic pairs increases on cooling below T l owing to chemical ordering.

It is important to note that the changes in AP1/AP2 ratio (APR) for the Zr 55 Cu 30 Ni 5 Al 10 liquid alloy from T l to T g (the ratio in Fig. 3 On cooling there is a change in partial pair distribution functions position fitted by the Gaussian function, background subtracted. The results for two major atomic contributions are shown in Table 1.

One can see the shift of Zr-Cu and Zr-Zr peaks to a higher R value on cooling. Intensification of both Zr-Zr and Zr-Cu peaks (area) is also observed. The Voronoi polyhedral analysis was also applied to the cells. However, owing to small cell sizes no systematic variation in coordination numbers of the polyhedral clusters is found with temperature except for Ni with CN=12. The number of Ni-centered clusters with CN=12 increased from 5 to 10 and to 12 on cooling from 1400 to 500 K. This is in line with intensification of the Ni PDOS peak in Fig. 6 indicating formation of Ni-centered clusters.

Conclusion

Fig. 1 .

 1 Fig. 1. (a) The experimental structure factor S(Q) at 1335 and 388 K calculated from the corresponding XRD intensity profile after the necessary corrections. (b) three typical G(R) functions and (c) three typical PDF(R) functions. Violet arrows (upward) show intensification of two sub-peaks responsible for Zr-(Cu,Ni) and Zr-Zr atomic interactions in the first coordination shell and two subpeaks at 0.5 and 0.58 nm in the second coordination shell. The arrow focusing downward-left indicates a decrease in the number of high R value interatomic distance on cooling.

Fig. 2 .

 2 Fig. 2. Fitting of four G(R) maxima (black curve) from 0.23 to 1.15 nm obtained at 501 K using seven Gaussian functions related to P1 -P7 (as indicated) (a). The baseline was subtracted prior to fitting according to the minima at about 2.2, 4, 6.7, 9 and 11.5 A. The inset shows R i 50 value at which two areas A1 and A2 (purple and green) are equal. The positions of P5 and P6 as a function of temperature (b). The intensity difference between the G(R) curves at 501 and 1335 K before the baseline correction (c).

Fig. 3 .

 3 Fig. 3. Positions of P1 and P2 in the first coordination shell (a) and the ratio of the corresponding peak areas (AP1/AP2) (b) as a function of temperature (a,b) and inverse temperature in (b). The positions of P3 and P4 as a function of temperature (c). The arrows in (a and c) indicate Yaxes corresponding to the data while the arrows in (b) indicate the corresponding X-axes. The insert in (a) is schematic representation of atomic redistribution within and between the first and higher order coordination shells leading to continuous structure changes in metallic liquids on heating and cooling leading to redistribution of the atomic number density as indicated by the red double side arrow.

  , in addition to intensification of 4 sub-peaks marked with arrows, there is a decrease in the fraction of high R value atomic pairs in the first coordination shell (from 0.35 to 0.4 nm, see downward-left arrow) on cooling which is missed when the baseline is corrected. Nevertheless the Gaussian fitting peaks positions represent the most probable interatomic distances. Also, although the positions of the maxima in the first coordination shell are strongly influenced by the asymmetry of the interatomic potentials (steeper slope at lower R) no shift is observed in the average atomic position ((P3+P4)/2) in the 2 nd coordination shell which is less affected by the potential. Also, the area under the first maximum (see Fig.2a) after the baseline correction was integrated and the point corresponding to 50 % of the integrated area (R i 50 ) (the areas of under the curve below (A1) and above this value (A2) are equal as shown in Fig. 2a,inset) was found in a wide temperature range below the liquidus temperature. As a result on cooling the R i 50 position changes insignificantly from 0.303 to 0.304 nm. The change in R i 50 is smaller than those obtained by the Gaussian fitting (Fig. 3a) but the tendency is same -larger R i 50 values are found at lower temperatures. On the other hand integration of the total area under P1 and P2 in the first coordination shell of the radial distribution function (RDF) at 501 K indicates the coordination number (CN) of 13.3, a local RDF maximum (RDF m ) at 0.317 nm and R i 50 =0.312 nm, while at 1335 K CN=13.0, RDF m = 0.313 nm and R i 50 =0.315 nm. Opposite variation in RDF m and R i 50 also indicate structural modification in the liquid phase on cooling. Owing to chemical and size similarity between Cu and Ni this quaternary alloy can be considered as a pseudo ternary Zr-(Cu,Ni)-Al one. Then, there are 6 atomic pairs in the studied alloy with distinct interatomic distances. The interatomic distances in the first coordination shell correspond to the sums of the Goldschmidt atomic radii of r Zr =0.160 nm, r Cu =0.128 nm, r Ni =0.125 nm and r Al =0.143 nm [84]. The first sub-peak (P1) in the first coordination shell corresponds to the nearest Zr-(Cu,Ni) distances (r Zr+Cu =0.288 nm) while the second sub-peak (P2) mostly corresponds to that of Zr-Zr pair (r Zr+Zr =0.320 nm) (Figs. 2a and 3a). The experimental peak position values are slightly lower than those predicted by the sum of the Goldschmidt atomic radii owing to formation of chemical bonds. Zr-Al interatomic distances are also closer to the value for the second peak and modify its position. These findings imply the possible structural changes and atomic redistributions within and between the coordination shells as schematically shown in the inset of Fig. 3a. Several clusters are located around the central one. The atoms which are considered as the solute-centering are black while surrounding atoms are white. Red circles schematically illustrate the first, second and third atomic coordination shells. Red arrow indicates permanent atomic number density redistribution between the corresponding shells on cooling. One should keep in mind that such a planar representation of the atomic structure is oversimplification and can be redrawn by using any other atom as the central one because any atom in the structure is a center of its own atomic cluster with a unique atomic arrangement, interatomic distances and a certain number of nearest neighbors. G(R) indicates statistically average atomic structure where maxima correspond to the most popular interatomic distances. Nevertheless this picture gives a visual representation of the structural processes occurring in a metallic fluid on cooling.

Fig. 4 .

 4 Fig. 4. G(R) calculated by using ab-initio MD simulation. The curves were smoothed over 10 points.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Partial PDF(R) functions obtained by computer simulation at 1400 K (liquid state) (a,b) and 500 K (glass state) (c,d).

  According to the results, strong chemical ordering forming Zr-Cu,Ni, Zr-Al and Zr-Zr atomic pairs takes place in the Zr 55 Cu 30 Ni 5 Al 10 supercooled liquid alloy on cooling below the liquidus temperature. However, here the change in the Zr-Cu,Ni peak area to other peaks area ratio (APR and APR n ) is smaller than in case of the Pd-Cu-Ni-P alloy studied earlier (Cu,Ni-P to other peaks) in accordance with a lower fragility index of the Zr 55 Cu 30 Ni 5 Al 10 melt. It is concluded that fragility is a sign of instability of short and medium range order in fragile liquids. MD simulation was found to support the observed structure variation as a function of temperature and indicated intensification of Zr-Cu, Zr-Al and Zr-Zr interactions on cooling.

  is 0.30/0.12=2.5) are smaller in the absolute values compared to those found in the Pd 40 Cu 30 Ni 10 P 20 liquid alloy (0.09/0.02=4.5)[START_REF] Louzguine-Luzgin | Structural basis for supercooled liquid fragility established by synchrotronradiation method and computer simulation[END_REF]. The AP1/AP2 ratio values normalized per Kelvin (APR n ) according to the supercooled liquid region on cooling (T l -T g ) are 0.005 (APR n =2.5/530 K -1 ) and 0.017 (APR n =4.5/270 K -1 ), respectively. This fact is in line with a lower fragility of the Zr 55 Cu 30 Ni 5 Al 10 melt compared to that of the Pd 40 Cu 30 Ni 10 P 20 one. Continuous structural changes in the supercooled liquid leading to the formation of atomic clusters likely change the activation energy for viscous flow, and thus, determine fragility of this glass-forming liquid through the equations (1 and 4).Fragility of liquids also correlates with the difference in the specific heat capacity (C p ) at T g between the liquid and glassy phases[START_REF] Klein | Excess thermodynamic properties of glassforming liquids: The rational scaling of heat capacities, and the thermodynamic fragility dilemma resolved[END_REF]. The Zr 55 Cu 30 Ni 5 Al 10 has C p values for liquid (C p ) phases at T g of about 38 and 28 J/molK, respectively, [ 86 ] while these values for the Pd 42.5 Cu 30 Ni 7.5 P 20 alloy are 45 and 28 J/molK, respectively [87]. Thus, C p l-g values of 10 and 17 J/molK for Zr 55 Cu 30 Ni 5 Al 10 and Pd 42.5 Cu 30 Ni 7.5 P 20 alloys, respectively, also correlate with a larger fragility of the latter alloy. As melting enthalpy of Zr 55 Cu 30 Ni 5 Al 10 is close to 7.5 kJ/mol and the liquidus temperature is about 1150 K [72,73], its melting entropy (ΔS m ) is about 6.5 J/Kmol. The resulted thermodynamic fragility m t parameter calculated using the formula: For an alloy close in composition to the Pd 42.5 Cu 30 Ni 7.5 P 20 one melting enthalpy is close to 5 kJ/mol [88] which gives ΔS m of 6.1 J/Kmol. The resulted m parameter is 104. If one uses ΔS m of 8.03 J/Kmol [89] then m=79 which is still higher than m=58 found for the Zr 55 Cu 30 Ni 5 Al 10 alloy.

		l ) and glassy
	(C p g
	m t =37.5ΔC p /ΔS m	(5)
	[85] is 58.

Table 1 .

 1 Position and area of two Zr-Zr and Zr-Cu peaks at two temperatures.

	Temperature, K		Zr-Zr		Zr-Cu
		R, nm	Area	R, nm	Area
	1400 K	0.3185	0.788	0.2850	1.022
	500 K	0.3225	0.850	0.2860	1.094

Acknowledgement:

This paper is dedicated to the memory of our wonderful colleague Alain Resa Yavari who started this research but passed away before it was completed. D.L. gratefully acknowledges financial support received from the World Premier International Research Center Initiative (WPI), MEXT, Japan. R.V.B. is grateful to the crew of Center for Computational Materials Science and E-IMR center at the Institute for Materials Research, Tohoku University, for continuous support.