Structure-Dynamics Relationship in Al-Mg-Si Liquid Alloys

Alaa Fahs and Philippe Jarry

Univ. Grenoble Alpes, CNRS, Grenoble INP, SIMaP, F-38000 Grenoble, France and
C-TEC, Parc Economique Centr'alp, 725 rue Aristide Bergès, CS10027, Voreppe 38341 cedex, France

Noël Jakse

Univ. Grenoble Alpes, CNRS, Grenoble INP, SIMaP, F-38000 Grenoble, France

(Dated:)

Enhancing properties and performances of aluminium alloys by a control of their solidification is pivotal in automotive and aerospace industries. The fundamental role of the structure-diffusion relationship is investigated for Al-Mg-Si liquid alloys taken as a prototype of Al-6xxx. For this purpose, first principles-based molecular dynamics simulations were performed for various Si and Mg content for Al-rich compositions, including the binary alloy counterparts. Results indicate that Mg and/or Si in alloys create a more compact ordering around Al in the alloy than the one existing in pure Al, thus lowering diffusion. Mg promotes icosahedral short-range order, while Si displays a preference towards cubic local ordering, impacting diffusion based on their respective content. It suggests a mechanism whereby an increase in Mg content generally lowers the diffusion of each species, whereas an increase in Si content enhances their diffusion, providing insights for future alloy design.

# I. INTRODUCTION

Aluminium alloys represent one of the main categories of structural metallic materials widely utilized in automotive construction and the aerospace industries [1–5]. They are attractive due to their high strength-to-density ratio, functional extrudability, age hardening characteristics, excellent corrosion resistance, as well as surface and welding properties. Moreover, the possibility of an almost complete recycling is a strong argument from an environmental point of view. Controlling properties and improving quality and performance of Al-based alloys represent timely and major industrial challenges [6, 7].

Improvements can be possible through controlled modification of the targeted microstructure. This requires a good understanding of solidification, a non-equilibrium process governed by the trade-off between solute diffusion in the liquid phase, crystal-liquid interfacial energy and anisotropy [8]. Unfortunately, experimental data for self-diffusion in liquid metals are scarce and most of them related to monatomic liquid metal, mainly due to a lack of specific radio-isotopes [9]. Alternatively, inserting the experimental viscosity data, widely available in the literature, into phenomenological laws such as the Stokes-Einstein (SE) equation [10, 11], could indirectly give a rough estimation of solute diffusivity. However, this approximated method happens to fail in most of the cases for multi-component liquid alloys. It might be due to a possible different local ordering around each species [12, 13], which falls out of the assumption of

the SE relation, and does not extend to undercooling conditions where solidification occurs.

Similarities in densities and coordination numbers between the liquid and crystal, for instance, led to the recognition that the local atomic ordering, known as the short-range order (SRO), might be different from the one leading to the long-range crystalline structure [14]. Frank [15] was the first to show that the icosahedral ordering is locally more stable. This might engender complex situations in the undercooled melt in terms of the variety of polymorphs [16], competing short range orders [17–19], or an interplay between chemical and fivefold symmetry (FFS) orderings in liquid alloys [20–22] prior to nucleation. More precisely, the proportion of atom pairs belonging to Icosahedral Short-Range Order (ISRO) strongly depends upon the nature of the solutes. In particular, in the Al₉₃Cr₇ liquid alloy, *ab initio* molecular dynamics (AIMD) calculations showed the existence of an onset of Icosahedral Medium-Range Order (IMRO), induced by the interplay between Chemical Short-Range Order (CSRO) and ISRO. This leads to dynamical heterogeneities (DHs) characterized by a strong heterogeneity in the propension to diffuse [21, 23]. Another AIMD study of the liquid Al₁-xNix alloy pointed out the existence of a non-linear evolution with composition x of dynamic properties due to non-linear composition dependence of CSRO [24]. More generally, a deeper understanding of the structure-dynamic relationship in liquid metals and alloys is needed for designing alloys with desired properties.

In this context, the possibility of calculating dy-
namic properties together with the underlying struc-
ture in the liquid and undercooled states on the same
footing can be provided by *ab initio* molecular dy-
namics [12, 13, 20, 21, 23–33]. For most metallic
systems, AIMD is able to accurately account for in-
teratomic interactions on the basis of the Density
Functional Theory (DFT) [34]. Despite the heavi-
ness of the underlying quantum computation to get
the electronic structure, AIMD simulations currently
allow us to build phase space trajectories of system
sizes of the order of 10² atoms over several hundreds
of picoseconds [21], relevant for the determination of
diffusion coefficients for all species. Another strat-
egy consists in developing Machine Learning (ML)
interatomic potentials trained on an appropriate set
of AIMD configurations, allowing to perform up to
million-atoms molecular dynamics (MD) simulations
with an accuracy close to *ab initio* [35]. It should
be noted that ML potentials for aluminium and its
alloys in the liquid and solid states were proposed
recently [36–38] and used to study early stage of so-
lidification [37].

The present work aims at investigating the
structure-dynamic relationships in liquid and un-
dercooled Al-Mg-Si ternary alloys, which are proto-
typical to the Al-based 6xxx used in the automo-
bile industry. For this purpose, *ab initio* molecular
dynamics are performed for various Al-rich compo-
sitions, namely, Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub>, Al<sub>70</sub>Mg<sub>20</sub>Si<sub>10</sub>, and
Al<sub>70</sub>Mg<sub>10</sub>Si<sub>20</sub>. In particular, the addition of Mg and
Si in liquid aluminium is further quantified by per-
forming also AIMD simulations for Al<sub>90</sub>Mg<sub>10</sub> and
the Al<sub>90</sub>Si<sub>10</sub> binary counterparts to disentangle ef-
fects of Mg and Si in the ternary liquid alloys. It
should be mentioned that these compositions have
a higher Si and Mg content than the 6xxx Alloys.
This choice was made given the small size of the sim-
ulation boxes, and intends to highlight the generic
relationship between local atomic ordering and the
dynamics. Besides, due to eutectic solute partition-
ing upon crystallization, the local liquid composition
is significantly higher at the solid-liquid interface
than in the liquid bulk, typically by a factor 3 for Mg
and 10 for Si, which makes the compositions studied
here relevant. Our results show significant structural
heterogeneities with Mg or Si addition in aluminium
either in binary or ternary alloys. In all cases, an
enhancement of the FFS around Mg atoms is ob-
served, while a decrease is seen around Si, leading to
a decoupling of diffusion. Impact of these findings
on the early stages of solidification in the alloys is
further discussed.

The layout of the paper is the following: Sec. II is
devoted to the computational aspects of the work, Sec.

III shows results of structural and dynamic properties
along with a discussion of the interplay between the
local structural properties and the dynamics, and
finally, Sec IV draws the conclusions.

II. COMPUTATIONAL BACKGROUND

A. Ab initio molecular dynamics

The structure and dynamic properties of bi-
nary Al<sub>90</sub>Mg<sub>10</sub> and Al<sub>90</sub>Si<sub>10</sub>, as well as ternary
Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub>, Al<sub>70</sub>Mg<sub>20</sub>Si<sub>10</sub>, and Al<sub>70</sub>Mg<sub>10</sub>Si<sub>20</sub> liq-
uid alloys, were determined by means of *ab initio*
molecular dynamics simulations within the density
functional theory, using the Vienna *Ab initio* Simu-
lation Package (VASP) [39]. Projected augmented-
wave method was used to describe the electron-ion
interaction [40], with a plane wave cutoff of 250 eV
and only the Γ-point to sample the Brillouin zone.
Exchange and correlation effects were taken into
account through the Local Density Approximation
(LDA)[41]. As a matter of fact, it was shown in pre-
ceding works that the LDA, among other functionals,
gives the best representation of the structure and dif-
fusion coefficients with respect to experiments in pure
Al and Al-rich alloys in the liquid and undercooled
states [23, 33, 42, 43].

The phase space trajectories were constructed
through dynamical simulations by numerical inte-
gration of Newton's equation of motion using the
standard Verlet algorithm [44] in the velocity form
with a time step of 1.5 fs. All simulations were done
in a cubic box of N = 256 atoms with periodic bound-
ary conditions in the three directions of space. The
number of each species in all the simulations can
be found in Table I of the Supporting Information
(SI) file [45] (see also references [33, 37, 43, 46-48]
therein). For each simulation, the volume V was
adjusted to reach pressure with an uncertainty less
than 0.5 GPa, which is of the order of pressure fluc-
tuations (see Fig. S1 in the SI file [45]). In order
to control the temperature T, the simulations were
performed in the canonical ensemble (constant N, V,
and T) with a Nosé thermostat [49]. Thermalization
was carried out for at least 20 ps, according to the
relaxation time evaluated from the self-intermediate
scattering function F<sub>S</sub>(q, t) [50], that are given in Ta-
ble 4 in the SI file [45]. Simulations were continued
during 60 to 160 ps with decreasing temperature, in
order to extract the structural and dynamic prop-
erties with a sufficient statistics. The temperature
evolution in the undercooled states was obtained by
quenching the system step-wise down to 600 K with a temperature step of 100 K from an equilibrated state at the previous temperature, resulting in an average quenching rate of 10¹³ K/s.

From equilibrated trajectories at each temperature, a regular sampling of 10 configurations was performed to obtain Inherent Structure (IS) configuration [51], from which a structural analysis is extracted by means of the Common-Neighbor Analysis (CNA) [52] using Python modules of the Open VIsualization Tool (OVITO) [53]. To this end, static relaxation simulations were performed for each of these configurations, performing an energy minimization procedure with the conjugate gradient method. This procedure is important to uncouple the vibrational motion from the underlying structural properties by bringing atoms to a local minimum in the potential-energy surface and reduce the noise in the CNA index calculations.

## B. Structural and dynamic properties

Structural properties were first investigated through the calculation of the total structure factor using the Debye equation for the purpose of comparing the simulation results to existing experimental measurements:

$$S(q) = \frac{1}{N} \sum_{j,k} b_j b_k \langle e^{i\mathbf{q}[\mathbf{r}_j - \mathbf{r}_k]} \rangle; \quad (1)$$

where $\mathbf{q}$ is the wave vector, $b_j$ and $\mathbf{r}_j$ represent respectively the neutron or X-ray scattering length, and the position of the atom $j$, and the angle brackets denote averaging over configurations sampled along the phase-space trajectory. The partial structural factors were calculated in the same manner within the Faber-Ziman formalism [46].

The partial pair-correlation functions $g_{ij}(r)$ were also considered to give a view of the local structural properties in the real space. These functions are written as:

$$g_{ij}(r) = \frac{V}{N_i 4\pi r^2} \lim_{\Delta r \to 0} \frac{n_{ij}(r + \Delta r)}{\Delta r}, \quad (2)$$

where $i$ and $j$ refer to two different types of atoms, $N_i$ is the number of atoms of type $i$, $n_{ij}(r)$ is the mean number of atoms of type $j$ around an atom of type $i$ in a spherical shell of radius $r$ and thickness $\Delta r$ centered on an atom $i$ taken as the origin. The

corresponding partial coordination numbers or the number of first nearest neighbors of a reference atom was calculated using the following equation:

$$N_{ij} = \rho_j \int_0^{r_{min}} 4\pi r^2 g_{ij}(r) dr, \quad (3)$$

where $\rho_j = N_j/V$ is the partial number density, and $r_{min}$ corresponds to the distance at the first minimum of the corresponding partial pair-correlation function $g_{ij}(r)$.

The local ordering surrounding each pair of atoms is further investigated through common-neighbor analysis (CNA), which is able to discriminate between various local topologies like cubic, namely fcc and bcc, hexagonal such as hcp, icosahedral ordering referred hereafter as FFS, or ISRO, as well as more complex environments like the Frank-Kasper polytetrahedral arrangements. The classification of each pair of atoms according to the number and topology of their common neighbors is done by means of a set of three indices in the Faken and Jonsson [52] nomenclature. For example, a bulk fcc crystal signature is 421, the hcp one contains the same pair as well as an additional 422 pair in equal proportion. The 444 pair together with the 666 one is a signature typical of the bulk bcc crystals with proportion 6/14 and 8/14, respectively. On the other hand, the 555 pair, corresponding to a pentagonal bi-pyramid, is characteristic of icosahedral order, with distortion thereof signaled by the presence of 544 and 433 pairs.

The diffusion of an atom of type $i$ can be tracked by means of the mean-square-displacement as function of time $t$, which is written as follows:

$$R_i^2(t) = \frac{1}{N_i} \sum_{k=1}^{N_i} \langle [\mathbf{r}_k(t+t_0) - \mathbf{r}_k(t_0)]^2 \rangle_{t_0}, \quad (4)$$

where the angular brackets represent an averaging over time origins $t_0$. Then, partial diffusion coefficients $D_i$ could be deduced from $R_i^2(t)$ as follows:

$$D_i = \lim_{t \to +\infty} \frac{R_i^2(t)}{6t}. \quad (5)$$

The diffusion coefficient were systematically corrected using the scheme of Yeh and Hummer [54] to take finite size effects into account, with an estimation of the viscosity of pure Al [33].

The individual atomic dynamics in liquids is also inspected through the self-intermediate scattering
function (SISF) and the non-Gaussian parameter (NGP). The former, which reads

$$F_S(\mathbf{q}, t) = \frac{1}{N} \left\langle \sum_{1=1}^{N} e^{-i\mathbf{q} \cdot [\mathbf{r}_1(t_0+t) - \mathbf{r}(t_0)]} \right\rangle, \quad (6)$$

provides insights into the temporal evolution of the spatial correlations of individual atoms and describes how the position of a particle correlates in time. $F_S(\mathbf{q}, t)$ starts at 1 when $t=0$ and decays to 0 at long time, atoms being totally uncorrelated with their initial positions. The value at $1/e$ is taken as a definition of structural relaxation time, $\tau_s$. The non-Gaussian parameter is defined as

$$\alpha_2(t) = \frac{3\langle R^4(t) \rangle}{5\langle R^2(t) \rangle^2} - 1, \quad (7)$$

where the angular brackets stand for the average over time origins of the total mean-square displacement, $R^2(t)$, and total mean-quartic displacement

$$R^4(t) = \frac{1}{N} \sum_{k=1}^{N} \left\langle \left[ \mathbf{r}_k(t+t_0) - \mathbf{r}_k(t_0) \right]^4 \right\rangle_{t_0} . \quad (8)$$

For a Gaussian diffusion process, i.e. homogeneous, this parameter is zero. Positive values indicate a deviation from Gaussian diffusive behavior and thus are interpreted as the presence of dynamic heterogeneities within the system.

### C. Comparison with experiments

First of all, the reliability of our *ab initio* calculations is assessed by comparing calculated total structure factors $S(q)$ with existing experimental data [47, 55]. Figs. 1 (a) and (b) show AIMD results for Al$_{90}$Mg$_{10}$ and Al$_{90}$Si$_{10}$ liquid alloys at temperatures $T = 1000$ K and $T = 900$ K, respectively. A good agreement is obtained with X-ray diffraction data [47, 55] regarding the position and amplitude of the first peak and subsequent oscillations. Moreover, Al-Mg and Mg-Mg partial structure factors in liquid Al$_{90}$Mg$_{10}$ are compared to experimental data of Waseda [47] in Fig. S2 in the SI file [45] within the Faber-Ziman formalism [46]. The Mg-Mg partial structure factor is well reproduced. This is all the more true for Al-Mg partial, especially regarding the peak positions, even if their amplitude is a bit underestimated. Such a result is particularly remarkable given the experimental difficulties to extract partials. To the best of our knowledge, no experimental data were reported till date in the literature for the ternary Al-Mg-Si liquid alloys.

### D. Comparison with machine learning potentials

Machine learning interatomic potentials were developed recently for the Al-Mg-Si ternary system [38] using a High Dimensional Neural Network Potential (HDNNP) [56]. Using the version named NN-19, extracted from the Materials Cloud [57], local structural and diffusion properties were carried out by molecular dynamics (MD) using LAMMPS in the same conditions as our AIMD calculations. The use of HDNNP was possible by means of the library-based implementation n2p2 by Singraber et al. [58].

Partial pair-correlation functions of Al$_{80}$Mg$_{10}$Si$_{10}$ liquid alloy from HDNNP at $T = 1000$ K are plotted in Figs. 2 (a) and (b). As a general trends, AIMD curves are well reproduced by the HDNNP, especially for Al-Al, Al-Mg and Al-Si pairs. Only the first peak of the Mg-Si partial shows an overestimated amplitude, and thus a higher Mg-Si affinity. More importantly, diffusion coefficients shown in Fig. S3 in the SI file [45], are in good agreement with AIMD in the liquid and undercooled states. This demonstrates the reliability of this potential and opens the way to larger scale simulations for early stage of solidification, with close to *ab initio* accuracy in the range of compositions studied here. The Al$_{90}$Mg$_{10}$ and Al$_{90}$Si$_{10}$ binary counterparts of the Al-Mg-Si ternary system were also investigated with the HDNNP. Fig. 2 (c) comparison with AIMD for Al$_{90}$Mg$_{10}$, where a good agreement is seen. However, the self-diffusion coefficients reported in Fig. S4 [45] show increasing departure with temperature, yet keeping the correct hierarchy $D_{Mg} < D_{Al}$. For Al$_{90}$Si$_{10}$ it was not possible to run stable MD simulations in the liquid state. The HDNNP for the ternary Al-Mg-Si ternary system is seemingly not completely transferable to binary counterparts. Training set should include the corresponding binary alloys, and the monoatomic metals if it were to be used for studying solidification phenomena.

## III. RESULTS

### A. Structural properties

The local structural properties are first considered through the partial pair correlation functions calculated from the AIMD simulations using Eq. 2. They are shown in Fig. 3 (a)-(f) for the Al$_{80}$Mg$_{10}$Si$_{10}$ ternary liquid alloy. The first peak position of each
<figure><img src="image_1.png" /><figcaption>FIG. 1. AIMD calculations of the total X-ray structure factors in liquid Al<sub>90</sub>Mg<sub>10</sub> (a) and Al<sub>90</sub>Si<sub>10</sub> (b) (solid lines), at temperature *T* = 1000 K and *T* = 900 K, respectively. Experimental data are shown with the symbols [47, 55] for comparison.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>FIG. 2. Partial pair-correlation functions of Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub> (a) and (b) as well as Al<sub>90</sub>Mg<sub>10</sub> (c) at *T* = 1000 K. HDNNP results (dashed lines) are compared to AIMD data (solid lines).</figcaption></figure>

partial reflects the first neighbor average distances,
and thus the bonding length between the correspond-
ing type of pairs. Going from the liquid state close
to the liquidus at T = 1000 K to the undercooled
states at T = 600 K, the first peak does not change
significantly in its position whatever the partial. Nev-
ertheless, its height increases visibly, showing that
the short-range order becomes more and more pro-
nounced and well defined. This is all the more true

for the second oscillation whose amplitude increases.
Moreover, for Al-Al, Al-Si and Al-Mg partials it splits
in two sub-peaks or displays a shoulder, a feature
that is often attributed to the reinforcement of the
FFS [21, 24, 32, 33, 59, 60]. Fig. 3(f) shows that
the amplitude of the g<sub>SiSi</sub>(r) remains very small and
close to unity, indicating that Si atoms are randomly
dispersed in the Al host liquid. Table II in the SI file
[45] gathers the first peak position value of all the
<figure><img src="image_3.png" /><figcaption>FIG. 3. Partial pair-correlation functions of Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub> liquid alloy for all studied temperatures by AIMD simulations. The green dashed vertical line indicates the first atomic distances of pairs in the Mg<sub>2</sub>Si cubic crystal (The Materials Project [48]: reference mp-1367).</figcaption></figure>

partials. It is worth mentioning first that addition of
 a small amount of Mg and Si in pure Al leaves the
 first nearest distance unchanged with r<sub>AlAl</sub> = 2.7 Å,
 close to pure liquid Al [33, 61, 62]. It also shows
 that r<sub>AIMg</sub> is intermediate between r<sub>AlAl</sub> and r<sub>MgMg</sub>.
 The height of the first peak of g<sub>AIMg</sub>(r) is also inter-
 mediate. The same holds for r<sub>AISi</sub>, but in that case,
 the first peak of g<sub>AISi</sub>(r) is higher than g<sub>AlAl</sub>(r) and
 g<sub>SISi</sub>(r), showing a slightly higher affinity between Al
 and Si.

In order to investigate a potential ternary effect,
AIMD simulations were performed for two binary
alloys, namely Al<sub>90</sub>Mg<sub>10</sub> and Al<sub>90</sub>Si<sub>10</sub> in the same
range of temperatures. The results for the partial
pair-correlation functions are shown in Figs. S7 and
S8, and the corresponding bond lengths are included
in Table II in the SI file [45]. The same trends on
the partials are seen, and only little change can be
observed with respect to the ternary alloy. Neverthe-
less, addition of Si in Al-Mg binary alloys leads to

the occurrence of a significant Mg-Si bonding in the
ternary alloy as can be seen in Fig. 3(e). A similarity
can be observed between the calculated Mg-Mg or
Mg-Si distances and those reported for the Mg₂Si
crystal structure, whose first neighbor distances are
indicated by the vertical dashed lines in Fig. 3 (d)-(f).
The bond-angle distribution in Fig. S9 of the SI file
[45] also suggests short-range ordering compatible
with precursors of this phase.

The resulting partial and total coordination num-
bers and their temperature evolution are represented
in Fig. 4 for Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub>. Only a modest increase
of the total coordination numbers around Al, Mg
and Si is observed with decreasing temperature be-
tween T = 1000 K and T = 600 K. This mirrors the
weak evolution of the first peak of the partial pair-
correlation function. Structural changes are better
revealed from the splitting of the partial-pair corre-
lation function at the level of the second peak. This
corresponds typically to a change of the orientational
<figure><img src="image_4.png" /><figcaption>FIG. 4. Total and partial coordination numbers as function of temperature in the ternary Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub> liquid alloy. Uncertainties in the values are of the order of 0.2.</figcaption></figure>

ordering as will be seen in the common-neighbor analysis below. The total coordination number of Al is: $N_c^{Al} = 12.4 \pm 0.2$, which is an intermediate value between $N_c^{Mg} = 15.0 \pm 0.2$ and $N_c^{Si} = 10 \pm 0.2$. Similar values for the total coordination numbers were found for the two binary liquid alloys (see Fig. S10 in the SI file [45]). For all liquid alloys (binary and ternary), the total coordination numbers show the following trend: $N_c^{Si} < N_c^{Al} < N_c^{Mg}$. Moreover, addition of 10% Mg and/or 10% Si leads to a significant increase of coordination number around Al with respect to pure Al, whose values of $N_c$ are between 11.6 and 11.8 in the present range of temperatures [33, 62]. This feature implies (i) an increase of the compacity $\eta_{Al} = \pi \rho r_{max,AlAl}^3/6$, where $r_{max,AlAl}$ the position of the first peak and $\rho$ is the atomic density as shown in Table II of the SI file [45], since the coordination sphere remains the same as mentioned above, and (ii) the possible appearance of a specific chemical short-range order. The following inequality, $N_c^{Si} < N_c^{Al} < N_c^{Mg}$, could indicate that, in general, the local structure around each species is different, characteristic of a structural heterogeneity. It is worth noting that Mg-Mg and Mg-Si partial coordination numbers are larger than 1, and for Al-Si it is equal to 1, whereas Si-Si coordination is almost negligible. This confirms that, in contrast with Si-Si bonding, Mg-Mg pairing is favored in both Al<sub>90</sub>Mg<sub>10</sub> and Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub>,

and that Si atoms exhibit preferential coordination with Mg, as expected by the analysis of the corresponding partial pair-correlation function. This is characteristic of the appearance of a CSRO, the latter being revealed by means of the Warren-Cowley [63] order parameter, $\alpha_i = 1 - N_c^{i-Al}/x_{Al}(N_c^{i-i} + N_c^{i-Al})$, giving $\alpha_{Si} = -0.06$ in Al<sub>90</sub>Si<sub>10</sub> and $\alpha_{Mg} = 0.04$ in Al<sub>90</sub>Mg<sub>10</sub>.

The SRO is further investigated by means of CNA [52] around each species for the binary and ternary liquid alloys. In Table I, results are gathered for two temperatures showing the most important bonded pairs determined from IS configurations. Five-fold symmetry is represented by the 555 pairs as well as the 554 ones, which are a distorted version of 555 [52]. A high degree of FFS characterizes the SRO of all alloys in the liquid at T = 1000 K, and becomes even more pronounced in the undercooled state at T = 600 K.

As a general observation, for all systems, the FFS around Mg is higher than around Al, while for Si it is lower. For Mg, the presence of 666 pairs could indicate the presence of bcc ordering when associated to 444 pairs. However, the proportion of 444 is too small, suggesting therefore Frank-Kasper polyhedra, mainly Z14 and Z15, to be more likely and consistent with a coordination number close to 15. For Si, the presence of 421 and 422 pairs is an indication of a significant cubic symmetry. The respective specific SRO around Mg and Si is reinforced in the undercooled regime. The 433 pairs can be associated either to very distorted 555 pairs or to distorted cubic symmetry. It generally decreases with decreasing temperature, indicating that the liquid becomes more structured. For both temperatures, the addition of Mg or Si atoms leads to a global increase of the FFS in each alloys by up to 10%, and sometimes even more around Al atoms with respect to pure liquid Al [33]. The increase of FFS around Al in alloys is consistent with a higher compacity, as mentioned above [15]. Table I also reports the presence of a small amount of 300 and 200 pairs in ternary alloys, also found in the Mg<sub>2</sub>Si phase (see Table S3 in the SI file [45]) and consistently with the typical distances found in the partial-pair correlation functions and bond-angle distributions. Regarding the solubility phase diagram, Mg<sub>2</sub>Si compound is in equilibrium with Al, and there is a quasi-binary line Al-Mg<sub>2</sub>Si at the Mg:Si ratio of 1.73. The two binary invariantectic reactions occur for 12.5 wt% Si (liquid $\rightarrow$ Al + Si) at 850 K, and for 34.5 wt% Mg (liquid $\rightarrow$ Al + Al<sub>8</sub>Mg<sub>5</sub>) at 723 K. The ternary invariant reactions are liquid $\rightarrow$ Al + Mg<sub>2</sub>Si + Si at 828 K at nearly 5 wt% Mg and 13 wt% Si on one end, and liquid
<table><caption>TABLE I. Common Neighbor Analysis for five Al-based liquid alloys calculated at T = 1000 K and T = 600 K. Results for pure liquid Al, reported in Ref. [33], are also shown for comparison. The typical error bar for the calculated pair abundances are of the order of 0.01.</caption><thead><tr><th>CNA (%)</th><th>Al [33]</th><th>Al<sub>90</sub>Mg<sub>10</sub></th><th>Al<sub>90</sub>Si<sub>10</sub></th><th>Al<sub>80</sub>Mg<sub>10</sub>Si<sub>10</sub></th><th>Al<sub>70</sub>Mg<sub>20</sub>Si<sub>10</sub></th><th>Al<sub>70</sub>Mg<sub>10</sub>Si<sub>20</sub></th></tr></thead><tbody><tr><td>T = 1000K</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[5 5 5] + [5 4 4]</td><td>26.9</td><td>38.1</td><td>45.1</td><td>29.3</td><td>12.6</td><td>32.9</td></tr><tr><td>[6 6 6]</td><td>1.40</td><td>2.83</td><td>5.39</td><td>2.70</td><td>0.53</td><td>2.33</td></tr><tr><td>[4 4 4]</td><td>1.40</td><td>2.24</td><td>1.03</td><td>2.04</td><td>2.05</td><td>1.82</td></tr><tr><td>[4 3 3]</td><td>20.0</td><td>25.0</td><td>21.1</td><td>23.6</td><td>25.9</td><td>24.5</td></tr><tr><td>[4 2 2]</td><td>13.9</td><td>11.2</td><td>9.58</td><td>13.3</td><td>13.6</td><td>13.0</td></tr><tr><td>[4 2 1]</td><td>8.00</td><td>6.34</td><td>5.71</td><td>6.05</td><td>5.74</td><td>5.34</td></tr><tr><td>[3 2 2]</td><td>4.10</td><td>3.34</td><td>2.52</td><td>4.36</td><td>6.79</td><td>3.90</td></tr><tr><td>[3 1 1]</td><td>10.6</td><td>6.84</td><td>4.86</td><td>10.3</td><td>19.2</td><td>9.48</td></tr><tr><td>[3 0 0]</td><td>-</td><td>1.44</td><td>1.10</td><td>3.00</td><td>5.66</td><td>2.69</td></tr><tr><td>[2 0 0]</td><td>-</td><td>0.66</td><td>0.62</td><td>1.52</td><td>4.50</td><td>1.23</td></tr><tr><td>T = 600K</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[5 5 5] + [5 4 4]</td><td>-</td><td>48.8</td><td>56.9</td><td>37.3</td><td>22.1</td><td>46.1</td></tr><tr><td>[6 6 6]</td><td>-</td><td>3.23</td><td>9.79</td><td>2.39</td><td>0.82</td><td>3.67</td></tr><tr><td>[4 4 4]</td><td>-</td><td>1.85</td><td>1.84</td><td>1.70</td><td>0.32</td><td>2.14</td></tr><tr><td>[4 3 3]</td><td>-</td><td>21.8</td><td>16.3</td><td>22.8</td><td>23.0</td><td>22.4</td></tr><tr><td>[4 2 2]</td><td>-</td><td>9.74</td><td>6.33</td><td>13.0</td><td>13.6</td><td>9.96</td></tr><tr><td>[4 2 1]</td><td>-</td><td>5.96</td><td>3.49</td><td>7.07</td><td>6.38</td><td>4.56</td></tr><tr><td>[3 2 2]</td><td>-</td><td>2.14</td><td>1.02</td><td>3.04</td><td>6.53</td><td>2.39</td></tr><tr><td>[3 1 1]</td><td>-</td><td>4.06</td><td>1.41</td><td>7.24</td><td>16.6</td><td>4.90</td></tr><tr><td>[3 0 0]</td><td>-</td><td>0.78</td><td>0.56</td><td>1.54</td><td>3.49</td><td>1.21</td></tr><tr><td>[2 0 0]</td><td>-</td><td>0.22</td><td>0.25</td><td>0.85</td><td>4.29</td><td>0.42</td></tr></tbody></table>

→ Al + Mg₂Si + Al₈Mg₅ at 722 K with 32.2 wt% Mg and 0.37 wt% Si on the other end. Note that the Al₁₀Mg₁₀Si alloy lies beyond the solubility limit, meaning that the primary phase during solidification would be the Mg₂Si intermetallic compound [64]. Nevertheless, it is difficult at this stage to claim the presence of precursors of the Mg₂Si phase that might form during slow solidification [65], given the small amount of these pairs present within the short time span of the simulations. This would require to perform large-scale MD simulations of the homogeneous nucleation process with a retrained ML potential [38], which is beyond the scope of the present paper.

Interestingly, addition of Mg in pure Al to form Al₉₀Mg₁₀ or Si to form Al₉₀Si₁₀ leads respectively to a global increase or decrease of ISRO. Then, adding Si in Al₉₀Mg₁₀ or adding Mg in Al₉₀Si₁₀ to form the ternary Al₈₀Mg₁₀Si₁₀ alloy follows the same trend, namely decreasing or increasing the FFS of the initial

corresponding binary alloy. Finally, in the ternary alloys, the same rule holds again, as the FFS is favored globally and around Al when the Mg content is increased at constant Si composition and disfavored when Si content is increased at constant Mg composition. Therefore, modifying the respective compositions of Mg and Si affects the local order in Al-Mg-Si alloy. More particularly, adapting ISRO might impact diffusion [21], thereby contributing to control solidification morphology.

## B. Structure-diffusion relationship

In order to determine self-diffusion coefficients for all the studied alloys, the mean-square-displacement $R_i^2(t)$ (Eq. 4) is calculated for each species $i$ on the production stage of the simulations after an equilibration of 30 ps. As can be easily seen on the log-log
plot of the $R_i^2(t)$ curves in Figs. S24-S28 in the SI file [45], this physical quantity is divided into two main parts, starting at small times by the ballistic motions and followed by a linear diffusive regime roughly after $t = 20$ ps for all alloys. In time spans between the two regimes, a slowing down of the MSD can be seen, corresponding to the cage effects, which reinforces with decreasing temperature. The self diffusion coefficient is obtained by a linear regression of $R_i^2(t)$ in the long-time diffusive part and by using Eq. 5.

Results of diffusion coefficients are drawn as a function of composition at $T = 1000$ K (see Fig. 5). Starting from pure liquid Al, which was studied in a previous contribution [33], binary liquid alloys are obtained by adding Mg or Si respectively for $Al_{90}Mg_{10}$ and $Al_{90}Si_{10}$, leading to a downshift of diffusion coefficients of Al for both alloys. From the structural analysis in the preceding section, the enhancement of the compacity around Al atoms together with the more pronounced ISRO might be at the origin of this slowing down for the two binary alloys. Nevertheless, Fig. 5 clearly shows that adding Mg or Si in liquid Al does not yield to the same consequences. ISRO is more pronounced in $Al_{90}Mg_{10}$ than in $Al_{90}Si_{10}$ leading to lower diffusion of Al atoms. Moreover, in $Al_{90}Mg_{10}$ the higher icosahedral ordering around Mg atoms leads to an even lower value of $D_{Mg}$. For $Al_{90}Si_{10}$, the local ordering is rather cubic and less compact around Si giving rise to higher Si diffusion. It is worth mentioning that from AIMD simulation studies, Si is known as a fast diffuser [66], giving rise to a more important decoupling in the binary alloys. For both binaries, a decoupling of the self-diffusion coefficients is observed, which is more pronounced in $Al_{90}Si_{10}$, and might be a consequence of structural heterogeneities, as revealed by the CNA analysis. Such a decoupling also reveals dynamic heterogeneities as will be seen below

The $Al_{80}Mg_{10}Si_{10}$ ternary alloy can be formed, as mentioned above, either by adding Si in $Al_{90}Mg_{10}$ or adding Mg in $Al_{90}Si_{10}$, in substitution of Al atoms. Fig. 5 shows that both paths are consistent with the mechanism unveiled above in which the increase of Mg content promotes ISRO and slows down the diffusion of each species, while the increase of Si composition induces cubic local ordering around Si and disfavors ISRO with the consequence of enhancing the diffusion of the species. This is all the more true for $Al_{70}Mg_{20}Si_{10}$ and $Al_{70}Mg_{10}Si_{20}$ obtained from $Al_{80}Mg_{10}Si_{10}$, respectively by increasing the Mg content at constant Si composition and increasing Si content and maintaining the Mg composition.

The global diffusion is marginally impacted by

<figure><img src="image_5.png" /><figcaption>FIG. 5. Diffusion coefficients for all the liquid alloys as function of composition at $T = 1000$ K. The open symbols represent the binary liquid alloys (Red: $Al_{90}Mg_{10}$; Blue: $Al_{90}Si_{10}$) and the closed symbols represent ternary liquid alloys (Green: $Al_{80}Mg_{10}Si_{10}$; Dark yellow: $Al_{70}Mg_{20}Si_{10}$; Magenta: $Al_{70}Mg_{10}Si_{20}$). Similar tendencies of diffusion coefficients have been found for other undercooling temperatures (see Figs. S29 and S30 in the SI file [45]). The cross represent the value for pure liquid Al from [33, 37]. The global diffusion coefficient is $0.50 \text{ Å}^2/\text{ps}$ for $Al_{80}Mg_{10}Si_{10}$, $0.45 \text{ Å}^2/\text{ps}$ for $Al_{70}Mg_{20}Si_{10}$, and $0.49 \text{ Å}^2/\text{ps}$ for $Al_{70}Mg_{10}Si_{20}$).</figcaption></figure>

the substitution of Al by Si and/or Mg atoms in the studied composition range, which is consistent with previous works [67], while the difference in self-diffusion of the individual species is significant. Such a difference reveals dynamic heterogeneities, which can be quantified using the non-Gaussian parameter, $\alpha_2(t)$. The curves of $\alpha_2(t)$ of all the alloys are drawn in Fig. 6(a) for the undercooled states at $T = 600$ K. They all display a significant positive amplitude, and thus reveal dynamic heterogeneities, as compared to those in the Inset determined at $T = 1000$ K, where only a small positive amplitude is seen, and mainly due to the cage effect that vanishes rapidly above 1 ps. The DHs then appear as a combined effect of the introduction of both Si and Mg elements. However, increasing the Mg content to 10% and 20% contribute the most to the DHs, with respect to $\alpha_2(t)$ for $Al_{90}Si_{10}$ displaying the smallest amplitude.

The corresponding self-intermediate scattering functions of all alloys are shown in Fig. 6(b), which give additional information on the structural relaxation process. At short times, typically below 0.1 ps, SISF shows a Gaussian-like behavior mirroring
<figure><img src="image_6.png" /><figcaption>FIG. 6. Non-Gaussian parameter for all the liquid alloys at T = 600 K (a). The Inset corresponds to the same curves determined at T = 1000 K. Self-intermediate scattering function for all the liquid alloys at T = 600 K (b). A wave number $q = 2.7Å^{-1}$ corresponding roughly to the position of the first peak of the structure factor $S(q)$ is considered. The Inset corresponds to $F_S(q, t)$ for Al$_{70}$Mg$_{20}$Si$_{10}$ along with the contributions of each species.</figcaption></figure>

the ballistic regime seen in the corresponding MSD. In the undercooled region at T = 600 K, SISF reveals two relaxation mechanisms. The short-time one, usually referred to as the so-called β-relaxation and characterized by a shoulder, corresponds to the cage effect in the MSD. The final exponential decay to zero from this shoulder corresponds to a slower relaxation process referred as the so-called α relaxation. It should be noted that above the liquidus at T = 1000 K β-relaxation is absent, which can be seen in Fig. S31 the SI file [45], as only this exponential decay exists, being an indication of a homogeneous diffusion regime. It appears clearly in the undercooled states that additions of Mg leads

to an increase of the structural relaxation time (α relaxation), whose values are given in the SI file [45]. It also gives rise to a more pronounced shoulder corresponding to an increase of the β-relaxation. The Inset of Fig. 6(b) illustrates the contribution of the individual species to the relaxation processes in the typical case of Al$_{70}$Mg$_{20}$Si$_{10}$, and shows that Mg and Al are at the origin of the β-relaxation, not seen for Si that relaxes much faster. This confirm the major role of Mg in the existence of dynamic heterogeneities in the ternary alloys studied here.

## IV. CONCLUSION

In conclusion, AIMD simulations were conducted in the liquid and undercooled states for Al-rich Al-Mg-Si ternary system, prototypical of 6xxx alloys, together with their binary counterparts, namely Al-Mg and Al-Si, with the goal of understanding the relationship between structure and diffusion. Our results show that the presence of Mg and/or Si leads to a local structure around Al significantly more compact in the alloys than in pure Al, which has the consequence of a global lowering of the diffusion. Moreover, an opposite influence of Mg or Si addition in the alloys on the liquid structure and dynamics was found. Whereas the presence of Mg in the alloys promotes ISRO, Si atoms show a tendency to disfavor ISRO and to form a cubic local ordering. These effects are shown to impact the diffusion of all species depending on the respective Mg and Si contents, and suggest the mechanism by which an increase of Mg content generally slows down the diffusion of each species, while the increase of Si composition enhances their diffusion. This demonstrates that chemical short-range order (CSRO) between Mg and Si is able to neutralize the ISRO-killing effect of Si atoms without annihilating the ISRO-promoting effect of Mg atoms. The evidence of such an interplay between topological and chemical short range orders provides significant insight into alloy designing principles, as far as solidification morphology selection is concerned. It may constitute the basis of any future investigation of the effect of other additional elements to the structural and dynamic properties of this system, with the possibility to fine tune the compositions to control diffusion mechanisms and thus the solidification process.
ACKNOWLEDGEMENT

We acknowledge the CINES and IDRIS un-
der Project No. INP2227/72914, as well as CI-
MENT/GRICAD for computational resources. This
work was performed within the framework of the
Centre of Excellence of Multifunctional Architec-
tured Materials CEMAM- ANR-10-LABX-44-01
funded by the "Investments for the Future" Pro-

gram. Funding by France-Relance under grant num-
ber 21-PRRD-0046-01 is also gratefully acknowl-
edged. This work has been partially supported by
MIAI@Grenoble Alpes (ANR-19-P3IA-0003) and in-
ternational SOLIMAT projet ANR-22-CE92-0079-
01. Discussions within the French collaborative
network in high-temperature thermodynamics GDR
CNRS3584 (TherMatHT) and in artificial intelligence
in materials science GDR CNRS 2123 (IAMAT) are
also acknowledged.

[1] R. Holmestad, R. Bjørge, F. J. H. Ehlers, M. Torsæter, C. D. Marioara, and S. J. Andersen, J. Phys. Conf. Ser. **371**, 012082 (2012).

[2] H. W. Zandbergen, S. J. Andersen, and J. Jansen, Science **277**, 1221 (1997).

[3] C. Ravi and C. Wolverton, Acta Mater. **52**, 4213 (2004).

[4] A. G. Frøseth, R. Høier, P. M. Derlet, S. J. Andersen, and C. D. Marioara, Phys. Rev. B **67**, 224106 (2003).

[5] P. Jarry and M. Rappaz, C. R. Phys. **19**, 672 (2018).

[6] T. P. D. Lassance, F. Scheyvaerts, Eng. Fract. Mech. **73**, 1009 (2006).

[7] M. Mageto, *Al-Mg-Si Alloys: Micostructure, Hardness and Tensile Properties* (LAP LAMBERT Academic Publishing, 2010).

[8] J. A. Dantzig and M. Rappaz, *Solidification* (EPFL Press English Imprint, 2016).

[9] T. Iida and R. Guthrie, *The Physical Properties of Liquid Metals* (Oxford science publications, Clarendon Press, 1988).

[10] W. Sutherland, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science **9**, 781 (1905).

[11] A. Einstein, Ann. Phys. (Berlin) **322**, 549 (1905).

[12] N. Jakse and A. Pasturel, Phys. Rev. B **94**, 224201 (2016).

[13] A. Pasturel and N. Jakse, Appl. Phys. Lett. **109**, 041904 (2016).

[14] D. Turnbull, J. Chem. Phys. **20** (1952).

[15] F. C. Frank, Proc. R. Soc. Lond. A **215**, 43-46 (1952).

[16] P. Ronceray and P. Harrowell, EPL **96**, 36005 (2011).

[17] N. Jakse and A. Pasturel, Phys. Rev. Lett. **91**, 195501 (2003).

[18] S. Becker, E. Devijver, R. Molinier, and N. Jakse, Phys. Rev. B **102**, 104205 (2020).

[19] U. R. Pedersen, I. Douglass, and P. Harrowell, J. Chem. Phys. **154**, 054503 (2021).

[20] N. Jakse and A. Pasturel, Phys. Rev. B **78**, 214204 (2008).

[21] A. Pasturel and N. Jakse, npj Comput. Mater. **3**, 33 (2017).

[22] C. Tang and P. Harrowell, J. Chem. Phys. **148**, 044509 (2018).

[23] A. Pasturel and N. Jakse, Appl. Phys. Lett. **110**, 121902 (2017).

[24] N. Jakse and A. Pasturel, J. Chem. Phys. **143**, 084504 (2015).

[25] G. Kurtuldu, P. Jarry, and M. Rappaz, Acta Mater. **61**, 7098-7108 (2013).

[26] G. Kurtuldu, Ph.D. thesis, Lausanne, EPFL (2014).

[27] G. Kurtuldu, P. Jessner, and M. Rappaz, J. Alloys Compd. **621**, 283-286 (2015).

[28] M. Rappaz, P. Jarry, G. Kurtuldu, and J. Zollinger, Metall. Mater. Trans. A **51**, 2651-2664 (2020).

[29] N. Jakse and A. Pasturel, J. Chem. Phys. **141**, 234504 (2014).

[30] N. Jakse and A. Pasturel, AIP Adv. **7**, 105212 (2017).

[31] M. Bouhadja and N. Jakse, J. Phys. Condens. Matter **32**, 104002 (2019).

[32] A. Pasturel, E. S. Tasci, M. H. F. Sluiter, and N. Jakse, Phys. Rev. B **81**, 140202 (2010).

[33] N. Jakse and A. Pasturel, Sci Rep **3**, 3135 (2013).

[34] J. Hafner, J. Comput. Chem. **29**, 2044 (2008).

[35] J. Behler, Chem. Rev. **121**, 10037 (2021).

[36] J. S. Smith, B. Nebgen, N. Mathew, J. Chen, N. Lubbers, L. Burakovsky, S. Tretiak, H. A. Nam, T. Germany, S. Fensin, et al., Nat. Commun. **12**, 1257 (2021).

[37] N. Jakse, J. Sandberg, L. F. Granz, A. Saliou, P. Jarry, E. Devijver, T. Voigtmann, J. Horbach, and A. Meyer, J. Phys. Condens. Matter **35**, 035402 (2022).

[38] A. C. P. Jain, D. Marchand, A. Glensk, M. Ceriotti, and W. A. Curtin, Phys. Rev. Mater. **5**, 053805 (2021).

[39] G. Kresse and J. Hafner, Phys. Rev. B **47**, 558 (1993).

[40] G. Kresse and D. Joubert, Phys. Rev. B **59**, 1758 (1999).

[41] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. **45**, 566 (1980).

[42] T. Bryk, T. Demchuk, N. Jakse, and J.-F. Wax, Front. Phys. **6**, 6 (2018).

[43] F. Demmel, L. Hennet, and N. Jakse, Sci. Rep. (Nature) **11**, 11815 (2021).

[44] L. Verlet, Phys. Rev. **159**, 98 (1967).
[45] See Supplementary Information file at [link] for some technical details; partial structural factors of Al₉₀Mg₁₀; comparison of structure-dynamics results obtained by AIMD and HDNNP; gαβ(r) of Al₉₀Mg₁₀ and Al₉₀Si₁₀; comparison of bond angle distributions in the liquid Al₈₀Mg₁₀Si₁₀ alloy and in the Mg₂Si solid crystal; Nc for the two binary liquid alloys; atomic distances and the corresponding Al compacity; CNA results obtained for different crystal phases and for all liquid alloys; total and partial relaxation time and SISF; R²(t) for all liquid alloys; temperature and composition dependence of D. The Supplementary Information file also contains Refs. [33,37,43,46–48].

[46] T. E. Faber and J. M. Ziman, Philos. Mag. **11**, 153 (1965).

[47] Y. Waseda, M. Ohtani, and K. Suzuki, Zeitschrift für Naturforschung A **28**, 1002 (1973).

[48] https://materialsproject.org/materials/mp-1367?chemsys=Mg-Si.

[49] S. Nosé, J. Chem. Phys. **81**, 511 (1984).

[50] K. Binder and W. Kob, *Glassy materials and disordered solids: An introduction to their statistical mechanics* (World scientific, 2011).

[51] F. H. Stillinger and T. A. Weber, Phys. Rev. A **25**, 978 (1982).

[52] D. Faken and H. Jónsson, Comput. Mat. Sci. **2**, 279 (1994).

[53] A. Stukowski, Model. Simul. Mater. Sci. Eng. **18**, 015012 (2009).

[54] I.-C. Yeh and G. Hummer, The Journal of Physical Chemistry B **108**, 15873 (2004).

[55] S. Mudry, A. Korolyshyn, and I. Shtablavyi, JPCS **98**, 012016 (2008).

[56] J. Fehler and M. Parrinello, Phys. Rev. Lett. **98**, 146401 (2007).

[57] A. C. P. Jain, D. Marchand, A. Glensk, M. Ceriotti, and W. A. Curtin, Materials Cloud Archive (2021).

[58] A. Singraber, T. Morawietz, J. Fehler, and C. Delago, J. Chem. Theory Comput. **15**, 3075 (2019).

[59] T. Schenk, D. Holland-Moritz, V. Simonet, R. Bellissent, and D. M. Herlach, Phys. Rev. Lett. **89**, 075507 (2002).

[60] C. P. Royall and S. R. Williams, Phys. Rep. **560**, 1 (2015).

[61] N. A. Mauro, J. C. Bendert, A. J. Vogt, J. M. Gewin, and K. F. Kelton, J. Chem. Phys. **135**, 044502 (2011).

[62] M. M. G. Alemany, L. J. Gallego, and D. J. González, Phys. Rev. B **70**, 134206 (2004).

[63] J. M. Cowley, Phys. Rev. **77**, 669 (1950).

[64] A. ALLOYS, Boston: Butterworths (1976).

[65] Y. Liu, S. Kang, and H. Kim, Mater. Lett. **41**, 267 (1999).

[66] N. Jakse and A. Pasturel, Phys. Rev. B **79**, 144206 (2009).

[67] J. Qin, X. Li, J. Wang, and S. Pan, AIP Advances **9** (2019).