# Fileset

[1-s2.0-S0022309323005793-main.pdf](https://mdr.nims.go.jp/filesets/d868f265-c463-40b6-ba61-3498921f55dd/download)

## Creator

Atsushi Tanaka, Atsuki Saito, Takashi Murata, [Ayako Nakata](https://orcid.org/0000-0002-3311-6283), [Tsuyoshi Miyazaki](https://orcid.org/0000-0003-3534-4404)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Large-scale DFT calculations of multi-component glass systems (SiO2)0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) : Accuracy of classical force fields](https://mdr.nims.go.jp/datasets/5091a091-4dd0-4373-9717-e05e19668a93)

## Fulltext

Large-scale DFT calculations of multi-component glass systems (SiO2)0.70(Al2O3)0.13( X O)0.17 ( X = Mg, Ca, Sr, Ba) : Accuracy of classical force fieldsJournal of Non-Crystalline Solids 625 (2024) 122714A0Contents lists available at ScienceDirectJournal of Non-Crystalline Solidsjournal homepage: www.elsevier.com/locate/jnoncrysolLarge-scale DFT calculations of multi-component glass systems(SiO2)0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) : Accuracy of classical forcefieldsAtsushi Tanaka a,∗, Atsuki Saito a, Takashi Murata a, Ayako Nakata b, Tsuyoshi Miyazaki ba Research and Development Group, Nippon Electric Glass Co., Ltd., 2-7-1 Seiran, Otsu, Shiga 520-8639, Japanb Research Center for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, JapanA R T I C L E I N F OKeywords:Multi-component glassFirst-principles calculationsLarge-scale DFTClassical force fieldsMolecular dynamics simulationsA B S T R A C TAlthough molecular dynamics (MD) simulation is a powerful tool for investigating the atomic-scale structuresof complex materials, several challenges limit their reliable and accurate application to multi-componentglass systems. The available force fields (FFs) that can treat many elements in a multi-component glass arelimited, and even if such a FF exists, its accuracy is suspicious due to the large variety and complexity ofchemical environments in these materials. First-principles calculations based on the density functional theory(DFT) are reliable, but prohibitively expensive with conventional methods. In this study, we use large-scaleDFT techniques and demonstrate that it is possible to perform efficient and accurate DFT calculations ofmulti-component glass systems, such as (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca, Sr, Ba), containing about1000–5000 atoms. From the results of large-scale DFT calculations, we evaluate the accuracy of some classicalFFs, and show that the accuracy for non-bridging oxygen atoms is very low especially when the Si–O distanceis short. Large differences in the distribution of Si–O–Si angles observed in the FF-MD and DFT-MD simulationsand the unique electronic structure in the case of 𝑋 = Mg are also discussed.1. IntroductionCommercial glass products have improved the quality of every-day life. Most products are made of multi-component glasses, andtheir compositions vary depending on their roles. For example, soda-lime glasses are used for windows because of transparency and highproduction efficiency. Some glass fibers are made of alkaline earthalumino-borosilicate glass for high insulations [1]. Cover glass sheetsfor electric devices are made of alkaline-containing silica glasses, whichcan be ion-exchanged to strengthen their surfaces.The advantage of mixing several components is the improvementof the characteristic properties of glass. However, the mechanisms toimprove these properties are still unclear and several possible mech-anisms have been proposed so far. The glasses are random-networkmaterials and it makes difficult to obtain the detailed atomic-scalestructures from experiments. Zachariasen [2] proposed a set of ruleson network structures of vitreous silica and complex glasses. Basedon Zachariasen’s model and experimental results of X-ray scattering,neutron scattering and NMR, the structures of vitreous silica and somebinary glasses have been clarified to some extent. For example, themedium-range structural order of vitreous silica was investigated [3],∗ Corresponding author.E-mail addresses: atstanaka@neg.co.jp (A. Tanaka), MIYAZAKI.Tsuyoshi@nims.go.jp (T. Miyazaki).and the effects of adding alkali-oxides on silica tetrahedral networkwere observed [4], which revealed that alkali-oxides break silicon–oxygen bonds and generate non-bridging-oxygen (NBO) atoms. As thenumber of NBO atoms increases, silica cannot maintain stable 3Dnetworks, and the viscosities are expected to decrease.These experimental approaches to glass structures are very usefulfor single or binary glass systems. Similar approaches were also appliedto multi-component glasses [5], but it was found that the structuresare more complex when other oxides, such as aluminum, boron, alkali,or alkali earth oxides are added. Because of the complexity and manyvariations of the coordination properties, experimental analyses of themulti-component glasses are challenging, compared with that single orbinary glass systems.To remedy the difficulties in the experimental analyses, computa-tional approaches are useful. Molecular dynamics (MD) simulationshave been used to predict structures and characteristic properties ofglasses. Some classical force-fields (FF) [6–8] are designed for multi-component glasses as well. However, the reliability of such classicalFF is often uncertain, since the experimental information is limited.It is expected that first-principles calculations based on the densityvailable online 25 November 2023022-3093/© 2023 The Authors. Published by Elsevier B.V. This is an open access ahttps://doi.org/10.1016/j.jnoncrysol.2023.122714Received 23 May 2023; Received in revised form 29 September 2023; Accepted 31rticle under the CC BY license (http://creativecommons.org/licenses/by/4.0/).October 2023https://www.elsevier.com/locate/jnoncrysolhttp://www.elsevier.com/locate/jnoncrysolmailto:atstanaka@neg.co.jpmailto:MIYAZAKI.Tsuyoshi@nims.go.jphttps://doi.org/10.1016/j.jnoncrysol.2023.122714https://doi.org/10.1016/j.jnoncrysol.2023.122714http://crossmark.crossref.org/dialog/?doi=10.1016/j.jnoncrysol.2023.122714&domain=pdfhttp://creativecommons.org/licenses/by/4.0/Journal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.wtpiip22obcd𝐻witDBAeftaoumtotga3ft𝜙wpAtofcBosu2nitinimla𝜙m𝑘ttl4functional theory (DFT) can solve this problem, because they pro-vide reliable and accurate information even for unknown materials.There are several attempts to develop reliable FFs, including the recentmachine-learning FFs, from the first-principles calculations of silicateglass or melt systems at various conditions [9–15]. Recently glassstructures simulated using DFT-MD were reported to be more reliablethan those using classical FFs [16,17]. Although DFT studies of multi-component glass systems are potentially useful, the difficulty is thatlarge systems containing many atoms are needed to reproduce theirregular network structure with medium range order. The compu-tational cost of DFT calculations is expensive, especially when thenumber of atoms is larger than 1000 [18]. Although there are alreadysome examples of large scale DFT calculations that clarify the complexstructures of glass systems [19,20], there are only a few reports of theDFT study of multi-component glasses [16,17,21,22], which generallyneed more atoms compared with single or binary glass systems.In this study, we use large-scale DFT techniques implemented inthe large-scale and linear-scaling DFT code CONQUEST to study theaccuracy of classical FFs for multi-component glasses. We demonstratethat reliable and accurate DFT calculations of multi-component glasssystems, (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca, Sr, Ba), containingabout 5000 atoms are possible using a large-scale DFT technique calledthe multisite method. A comparison of the structures from the DFT-MDand FF-MD simulations and the analysis of the electronic structure of(SiO2)0.70(Al2O3)0.13(𝑋O)0.17 are also provided.The rest of the paper is organized as follows: In Section 2, anoverview of large-scale DFT techniques and the accuracy of the basis setused in the DFT calculations are provided. The classical FFs employedin this study and the preparation method to generate the snapshotstructures for benchmark tests are also explained. In Section 3, theresults of benchmark tests and large-scale DFT calculations of multi-component glass systems are presented. A summary of the results isgiven in Section 4.2. Computational methods2.1. Classical force fieldsIn this study, three kinds of potentials are used for the classical FFs.The first was proposed by Pedone et al. [6], called FF-Pedone hereafter.Its potential is expressed as:𝑈𝑖,𝑗 (𝑟) = 𝐷𝑖𝑗 [{1 − 𝑒−𝑎𝑖𝑗 (𝑟−𝑟0)}2 − 1] +𝐶𝑖𝑗𝑟12+𝑧𝑖𝑧𝑗𝑒24𝜋𝜖0𝑟, (1)where, 𝑟 is the distance between atoms 𝑖 and 𝑗. The first, second, andthird terms represent the short-range Morse function, repulsive contri-bution, and long-range Coulomb potential, respectively. For compari-son, we also used two other FFs that carry different sets of Buckinghamtype parameters, as shown in Eq. (2); proposed by Deng et al. [7], andWang et al. [8], referred to as FF-Deng and FF-Wang, respectively.𝑈𝑖,𝑗 (𝑟) = 𝐴𝑖𝑗 exp(−𝑟𝜌𝑖𝑗) −𝐶𝑖𝑗𝑟6+𝑧𝑖𝑧𝑗𝑒24𝜋𝜖0𝑟, (2)here the first, second, and third terms represent the short-range elec-ronic repulsion, van der Waals interactions, and long-range Coulombotential. The Buckingham type potential is known to have drawbacksn the short range contact region; therefore, the first and second termsn Eq. (2) are replaced by other FFs. The parameters of these FFs arerovided in Supplementary..2. Large-scale DFT methods.2.1. Local basis sets and support functionsIn the large-scale DFT code CONQUEST[23–26], the Kohn–Shamrbitals for the state 𝑛, 𝜓𝑛(𝐫), or the density matrix 𝜌(𝐫, 𝐫′) are expressedy local functions around each atom 𝑖 with the orbital 𝛼, 𝜙 (𝐫),2𝑖𝛼alled support functions. The Kohn–Sham orbitals can be obtained byiagonalizing the Hamiltonian matrix,𝑖𝛼,𝑗𝛽 = ∫ d𝐫𝜙𝑖𝛼(𝐫)𝐻̂𝜙𝑗𝛽 (𝐫) (3)ith the overlap matrix, since the support functions are not orthogonaln general. In Eq. (3), 𝐻̂ is the Kohn–Sham Hamiltonian operator usinghe pseudopotential method. For the exchange–correlation term in theFT, the generalized gradient approximation (GGA) with the Perdew–urke–Ernzerhof (PBE) functional [27] was used throughout this study.lthough it was reported that dispersion corrections can improve thenergetics among the various polymorphs of SiO2[28], we have con-irmed that the effects of the empirical dispersion corrections [29] onhe atomic forces in the present glass systems were almost negligiblend the corrections were not considered in this study.When considerably large systems containing more than a few tensf thousands of atoms need to be calculated, a linear-scaling mode issed with the density matrix minimization method. In the linear-scalingode, the density matrix, 𝜌(𝐫, 𝐫′), is optimized to minimize the DFTotal energy without calculating the Kohn–Sham orbitals.The support functions themselves are represented in terms of onef two basis sets: the pseudo atomic orbitals (PAOs); and blip func-ions [30]. In this study, we use PAOs as a basis set, which areenerated by the MakeIonFiles code provided in the CONQUEST pack-ge [31], using the pseudopotentials generated by ONCVPSP code [32,3] with the parameters given in PseudoDojo[34,35]. The supportunctions of atom 𝑖 are usually expressed by a linear combination ofhe PAOs of the atom, as shown in Eq. (4):𝑖𝛼(𝐫) =∑𝜇𝑐𝑖𝛼(𝜇)𝜒𝑖𝜇(𝐫), (4)here 𝜒𝑖𝜇(𝐫) is the PAO of the orbital 𝜇 of the atom 𝑖. It is alsoossible to use a PAO for each support function (primitive basis set).lthough a systematic convergence is difficult with the PAO basis set,he accuracy of the calculations can be improved if the number ofrbitals is increased and the basis set of double-zeta plus polarizationunctions (DZP) is accurate enough in most cases [36–38].Table 1 shows the optimized lattice parameters and bulk modulusalculated with the DZP basis sets for SiO2, Al2O3, MgO, CaO, SrO andaO. The results show that the optimized structures are similar to thosebtained using other DFT calculations, confirming that the DZP basisets used herein are accurate enough for these oxides. Therefore, wesed the DZP basis sets throughout this study..2.2. Multisite method for large-scale DFT calculationsSince the computational cost (CPU time) to diagonalize the Hamilto-ian in Eq. (3) scales with the cubic of the number of support functions,t is essential to decrease the number of support functions to reducehe cost of the calculations for large-scale systems. Note that the costn the linear-scaling calculations is also proportional to the cube of theumber of support functions. If we are to express the support functionsn Eq. (4), we need at least one support function for each angularomentum [37]. However, if we construct the support functions as theinear-combination of PAOs on a target atom, as well as its neighboringtoms, defined by a cutoff radius 𝑅MS,𝑖𝛼(𝐫) =𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑢𝑟𝑠∑𝑘∑𝜇∈𝑘𝐶𝑖𝛼(𝑘𝜇)𝜒𝑘𝜇(𝐫). (5)Then, the number of support functions can be reduced to that of ainimal basis set [42,43]. Here, 𝜒𝑘𝜇 is a PAO of its neighboring atoms(including 𝑖 itself) with the orbital 𝜇 and the support functions arehus called multisite support functions (MSSFs).In CONQUEST, two methods have been implemented to determinehe linear-combination coefficients of MSSFs. One of the methods isocal-filter-diagonalisation (LFD) proposed by Rayson and Briddon [42,4,45], and the other is numerical optimization [43]. In the LFDJournal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.mH𝑅couTable 1Optimized lattice parameters and bulk modulus for various oxides calculated with the DZP basis set. Experimental values are also presentedfor comparison. SiO2 (𝛼-quartz) and 𝛼-Al2O3 have hexagonal symmetry, while the other four oxides have a cubic (NaCl) structure.SiO2(𝛼-quartz) Al2O3(𝛼-Al2O3) MgO CaO SrO BaOLattice parameter: 𝑎opt , 𝑐opt 4.90, 5.39 4.74, 13.13 4.28 4.85 5.21 5.61Lattice parameter: 𝑎exp, 𝑐exp 4.92, 5.41a 4.67, 12.95b 4.22c 4.82c 5.16c 5.52cBulk modulus : opt (GPa) 46.6 230.6 148.6 112.1 81.2 65.3Bulk modulus : exp (GPa) 38a 264b 160c 111c 89c 74ca Ref. [39].b Ref. [40].c Ref. [41].ethod, local molecular orbitals are first calculated by solving the localamiltonian of a cluster comprising of the atoms within the range ofLFD from a target atom. Subsequently, the linear-combination coeffi-ients are obtained by projecting the occupied local molecular orbitalsnto trial vectors which are localized around the target atom. Whensing the LFD method, the 𝑅MS should be equal or shorter than thesubspace region in the LFD method 𝑅LFD. On the other hand, the linear-combination coefficients are directly optimized to minimize the totalenergy in the numerical optimization method. When the numericaloptimization method is used, the LFD method is first employed toobtain the initial guess of the linear-combination coefficients. If both𝑅LFD and 𝑅MS are large, these two methods would result in almost thesame accuracy. Since the MSSFs depend on the charge density of thesystem, which in turn depends on the MSSFs, the linear-combinationcoefficients need to be determined self-consistently [42].The accuracy of the MSSFs depends on the 𝑅MS and the calculationswith MSSF can reproduce the accuracy of the primitive PAO basisset if 𝑅MS is large enough. One of the advantages is that the differ-ence between the MSSF and primitive basis set calculations decreasesexponentially not only in gapped systems, but also in metallic sys-tems. In Section 3.2, we investigate the accuracy of MSSF calculationsin the multi-component glass systems and show that accurate DFTcalculations of about 5000-atom systems are possible with the MSSFmethod.The reduction of the computational cost by the MSSF method issignificant in the calculations of large systems, for which the cost ofthe LFD becomes negligible compared with the cost of diagonalizing𝐻𝑖𝛼,𝑗𝛽 in Eq. (3). The number of DZP basis functions for Si or O is 13,while the number of their support functions in the MSSF method is 4.As a result, the reduction of the required memory size and cpu timeby the MSSF method is estimated to be about 10.5 (= (13∕4)2) and 34(= (13∕4)3) folds, respectively. For more accurate basis sets, such asTZTP, the reduction rate is even higher: 45 for memory requirementand 308 for the computational time.2.2.3. Modeling of the glass structures by melt-quench processThe glass structures are generated by the MD simulations of melt-quench processes using classical FFs with the LAMMPS package. For allMD simulations, the time step is fixed at 1.0 fs. The initial atomic po-sitions for FF-Pedone and FF-Deng are generated by random numbers,while the structure made by the melt-quench process with FF-Pedonein the following paragraph, is used as the initial structure for themodel with FF-Wang to prevent the so-called Buckingham catastrophe.The size of the simulation cells for the initial structures with FF-Pedone and FF-Deng are determined by the experimental densities. InSection 3.1, three systems are investigated for the benchmark test ofthe classical FFs; pure SiO2 (vitreous silica), (SiO2)0.5(Al2O3)0.5, and(SiO2)0.5(CaO)0.5. The experimental densities of the three systems arereported to be 2.20, 2.74 [46,47] and 2.90 [48] g/cm3, respectively.These values at room temperature are calculated using the Archimedesmethod.In the melt-quench process with FF-Pedone, the glass systems arefirst melted at 6000 K for 500 ps using the NVT ensemble. Subse-quently, they are cooled down to 300 K at a rate of 5 K/ps in the3NVT simulation. After cooling, they are relaxed for 200 ps in the NVTsimulations at 300 K. In the case of FF-Deng, the systems are firstequilibrated at 300 K for 60 ps, then melted at 6000 K for 100 ps, andequilibrated at 5000 K for 100 ps. Following the equilibration of theliquid state, the systems are cooled down to 300 K at a cooling rateof 5 K/ps. These processes were conducted using the NVT simulations.Subsequently, the systems were relaxed at 300 K for 60 ps in the NPTensemble at zero pressure, followed by the NVT simulations at 300 Kfor 60 ps. Regarding the process with FF-Wang, the systems are firstmelted at 3000 K for 10 ps in the NVT simulations, followed by the NPTsimulations at 3000 K and zero pressure for 100 ps. They are cooleddown to 300 K at a cooling rate of 1 K/ps in the NPT simulations atzero pressure. Finally, they are relaxed at 300 K for 100 ps in the NPTsimulations at zero pressure, then NVT simulations at 300 K for 100 ps.The densities of the glass structures generated with FF-Deng and FF-Wang are 2.29 and 2.47 g/cm3 for vitreous silica, 2.77 and 3.03 g/cm3for (SiO2)0.5(Al2O3)0.5, and 2.88 and 3.00 g/cm3 for (SiO2)0.5(CaO)0.5,respectively. The results agree with the experimental values for thestructures made using the FF-Deng process. The densities of the struc-tures made using the FF-Wang process are higher than the experimentalvalues by about 10% for the first two systems, but the densities agreemuch better for (SiO2)0.5(CaO)0.5. Note that only NVT ensemble simu-lations are used in the melt-quench processes with FF-Pedone, and theexperimental densities cannot be reproduced using NPT simulations.3. Results3.1. Benchmark tests for classical force fields in vitreous silica and binaryglassesIn this section, we investigate the accuracy of the three types of clas-sical FFs, FF-Pedone, FF-Deng and FF-Wang, explained in Section 2.1,for vitreous silica and binary glasses. Three systems were used forthe benchmark tests; pure SiO2 (vitreous silica), (SiO2)0.5(CaO)0.5, and(SiO2)0.5(Al2O3)0.5. The number of atoms is 999 for vitreous silica and1000 for the two binary glass systems. The snapshot structures wereprepared by the melt-quench processes, explained in Section 2.2.3. Self-consistent DFT calculations were performed using a randomly selectedsnapshot structure for each system and for each FF. The calculatedatomic forces by DFT, 𝐹DFT, were compared with those by classicalFFs, 𝐹FF, and the three components of the atomic forces on the Si andO atoms in the vitreous silica system are plotted in Fig. 1(a) and (b),respectively. Fig. 1(c) and (e) show the forces on O and Ca atoms,respectively, in (SiO2)0.5(CaO)0.5, while Fig. 1(d) shows the forces onAl atoms in (SiO2)0.5(Al2O3)0.5. The correlation coefficient, 𝑅, and themean absolute deviation (MAD) between the atomic forces calculatedby the three FFs and DFT are presented in Table 2.The tendencies of the atomic forces on Si, Al, and Ca are shownin Fig. 1(a), (d), and (e). The forces on these atoms calculated by FF-Pedone qualitatively agree with those by DFT, although the amplitudeof the forces by FF-Pedone is systematically underestimated, resultingin the high 𝑅 value with large MAD values. FF-Deng basically showsa similar behavior with FF-Pedone. The agreement is a little betterin vitreous SiO2, worse in (SiO2)0.5(Al2O3)0.5, and very similar inJournal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 1. Comparison of the components of atomic forces calculated by classical force fields 𝐹FF and DFT 𝐹DFT, for (a) Si and (b) O atoms in SiO2, (c) O atoms in (SiO2)0.5(CaO)0.5,(d) Al atoms in (SiO2)0.5(Al2O3)0.5, and (e) Ca atoms in (SiO2)0.5(CaO)0.5. Three types of classical force fields, FF-Pedone, FF-Deng and FF-Wang are investigated.Table 2Correlation coefficient 𝑅 and the mean absolute deviation (MAD) between the atomic forces calculated by classical force fields and DFT.FF-Pedone FF-Deng FF-Wang𝑅 MAD 𝑅 MAD 𝑅 MADSi: SiO2 0.798 0.913 0.836 0.657 0.780 0.680O: SiO2 0.675 0.717 0.712 0.657 0.771 0.471Si: (SiO2)0.5(Al2O3)0.5 0.866 0.791 0.730 0.956 0.648 1.105O: (SiO2)0.5(Al2O3)0.5 0.727 0.655 0.554 0.801 0.772 0.493Al: (SiO2)0.5(Al2O3)0.5 0.880 0.515 0.717 0.586 0.727 0.796Si:(SiO2)0.5(CaO)0.5 0.811 0.918 0.832 0.643 0.721 0.841O:(SiO2)0.5(CaO)0.5 0.464 1.295 0.435 1.211 0.610 0.669Ca:(SiO2)0.5(CaO)0.5 0.613 0.436 0.518 0.434 0.526 0.594(SiO2)0.5(CaO)0.5. Using FF-Wang, the agreement with DFT forces invitreous silica is worse than those of the other two FFs. In addition,agreements in Ca atoms by these FFs are worse than those for otherelements in common.For the atomic forces on oxygen atoms, these three FFs show ratherdifferent accuracies. As seen in Fig. 1(b) and (c), the agreement is verydifferent between vitreous silica and (SiO2)0.5(CaO)0.5 when using FF-Pedone or FF-Deng. The agreement in the case of (SiO2)0.5(CaO)0.5 ismuch worse, as can be seen in Table 2. On the other hand, with FF-Wang, such a large difference cannot be seen between oxygen atomsin vitreous silica and those in (SiO2)0.5(CaO)0.5 Although it is difficultto judge which FF is more accurate in general, the three classicalFFs evidently have different characteristics, especially for the oxygenatoms.Aiming to clarify the conditions where 𝐹FF deviates significantlyfrom 𝐹DFT for oxygen atoms, the comparison of the forces dependingon the distance is investigated. Fig. 2 shows |𝐹FF|, |𝐹DFT|, and |𝛥𝐹 |(𝛥𝐹 = 𝐹FF − 𝐹DFT) for the oxygen atoms in (SiO2)0.5(CaO)0.5. Here,the oxygen atoms are classified into bridging oxygen (BO) and NBOatoms. Usually, Ca atoms exist near the NBO atom, on the other side ofthe O–Si bond. In Fig. 2, the horizontal axis denotes the interatomicdistance between an oxygen atom and its nearest neighboring atom. ForFF-Pedone and FF-Deng, the results show that the deviations between4FF and DFT forces, |𝛥𝐹 |, for NBO atoms are larger than those for BOatoms, implying that the accuracy of the FFs for NBO atoms is muchworse compared with the one for BO atoms. The deviation is especiallylarge when the bond length of NBO atoms is shorter than 1.6 Å. Incontrast, such large differences between NBO and BO atoms are notobserved in the atomic forces calculated by FF-Wang. Furthermore,the agreement with the DFT forces is better when using FF-Wang, inthe case of (SiO2)0.5(CaO)0.5, compared with FF-Pedone and FF-Deng.These results suggest that FF-Pedone and FF-Deng tend to underes-timate the repulsive potential for the NBO atoms having short Si–Obonds. The distance between oxygen and silicon atoms is sometimesshorter than 1.6 Å in the trajectory of MD simulations with FF-Pedoneand FF-Deng, while such short O–Si bonds are rarely observed withFF-Wang. This probably comes from the difference in the propertiesof NBO atoms between the first two FFs and FF-Wang, in the case of(SiO2)0.5(CaO)0.5.Generally, the bonds between Al and O atoms are mainly covalent,and BOs are dominant in (SiO2)0.5(Al2O3)0.5 glass. This may be thereason the agreements in (SiO2)0.5(Al2O3)0.5 are better than those in(SiO ) (CaO) , using the first two FFs.2 0.5 0.5Journal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 2. (a) Comparison of the atomic forces calculated by DFT (green crosses) and classical force fields (purple open circles) for O atoms in (SiO2)0.5(CaO)0.5. (b) Difference inthe two forces for the bridging oxygen (BO) and non-bridging (NBO) oxygen atoms, expressed by blue open circles and red crosses, respectively. In all the figures, the horizontalaxis shows the shorter distance of oxygen atoms with the neighboring atoms.Fig. 3. Distribution of the number of Ca atoms having various coordination numbers in (SiO2)0.70(Al2O3)0.13(CaO)0.17 systems, containing (a) 216, (b) 927, (c) 4944, and (d) 9888atoms. The snapshot structure of the 4944-atom system is shown in (e).3.2. Large scale DFT calculations of the multicomponent glass systems,(SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca, Sr, Ba)Next, we analyzed the accuracy of classical FFs for more com-plicated, but widely used multi-component glass systems (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 where group II element 𝑋 is Mg, Ca, Sr, or Ba.Here, we used FF-Pedone because the parameters for various 𝑋s areavailable only with this FF. For such multi-component glasses, the sizeof the system should be large enough to properly reproduce the actualglass structures. Fig. 3(a)–(d) show the distribution of the coordinationnumber around Ca atoms in the MD simulations for different sizes of thesystem. This analysis suggests that we need about 5000 atoms to repro-duce the correct distribution of the coordination number of Ca atoms. Asnapshot structure of MD simulations of (SiO ) (Al O ) (CaO) ,52 0.70 2 3 0.13 0.17which contains about 5000 atoms, is shown in Fig. 3(e). If the atomicforces depend on the distribution of the coordination numbers, it isimportant to calculate the electronic structure and atomic forces ofsuch large systems by DFT. As discussed in 2.2, it is very expensive toperform accurate DFT calculations of such large systems using the con-ventional DFT methods; therefore, we demonstrate here that efficientand reliable calculations are possible with the MSSF method.Fig. 4 shows the accuracy of the MSSF method in the calculation ofvitreous silica glass systems. As explained in Section 2.2.2, the accuracyof the method is controlled by the multisite range, 𝑅MS. The totalenergy of the vitreous silica system containing 999 atoms calculated bythe MSSF method using different values of 𝑅MS is presented in Fig. 4(a).For this size, the exact value can be evaluated using the primitive PAObasis sets, which is shown by the dotted red line in the figure. TheseJournal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 4. (a) Total energy of a vitreous silica system containing 999 atoms, calculated by the MSSF method using different multisite range 𝑅MS. Comparison of the calculated atomicforces of the snapshot structure using the primitive PAO basis set and the MSSF calculations using (b) 𝑅MS = 4.2 Å and (c) 𝑅MS = 8.5 Å.results show that the MSSF calculation with 𝑅MS = 6.4 Å is already veryaccurate. The accuracy of the forces calculated using the MSSF methodwith 𝑅MS = 4.2 and 8.5 Å is also investigated, and the results are shownin Fig. 4(b) and (c), respectively. The agreement of the calculated forcesby the MSSF method with those using a primitive basis set is almostperfect for both cases.Next, the DFT calculations of the systems containing about 5000atoms (Fig. 3(e)) using the MSSF method with 𝑅MS = 8.5 Å arepresented. The results in Fig. 4 show that the calculation with thisrange of 𝑅MS must be highly accurate. It was confirmed that the SCFcalculations of such large systems with the MSSF method are stableand robust. The calculations were performed using 9–54 nodes (1 node= 2 CPUs of Xeon Platinum 2.9 GHz, 24 cores) of the NumericalMaterials Simulator at NIMS. The elapsed time to converge the SCFcharge density and calculate the atomic forces was 6367 s using 9nodes, and 2298 s using 36 nodes. It is still expensive to perform long-time MD simulations with this timing. However, it should be noted thatthe elapsed time can be significantly reduced by using 𝑅MS = 4.2 Å.Furthermore, the number of SCF iterations is expected to be muchsmaller after the first MD step, since we can begin with the MSSFcoefficients and self-consistent charge density obtained in the last MDstep. Such an approach will be studied in the future.The calculated atomic forces are compared with those by FF-Pedonein Fig. 5, showing the |𝛥𝐹 | for the BO and NBO atoms. These BO andNBO atoms are further classified by their nearest neighbor atoms, Sior Al. Again, the agreement of the forces between DFT and FF-Pedoneare poorer for the NBO atoms than the BO atoms. Furthermore, for theoxygen atoms bonded with Si atoms, the difference in the forces is verylarge when the distance is shorter than 1.55 Å. This behavior is commonfor all cases of atom 𝑋.As earlier mentioned, various kinds of local structures can be re-produced in the large-scale structural models containing about 5000atoms. It is expected that more revised FFs, including the machine-learning FFs, would be developed and applied to such large-scalemulti-component glass systems in the future. Here, it is demonstratedthat accurate DFT calculations can be performed for the complex glasssystems containing more than 5000 atoms, and the accuracy of futureFFs can be clarified in detail by such large-scale DFT calculations.3.3. DFT-MD simulations of (SiO2 )0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 =Mg, Ca,Sr, Ba)Next, we discuss the MD simulations based on DFT. Although it ispossible to perform MD simulations of 5000-atom systems with theMSSF method in principle, it is still expensive. Thus, the DFT-MDsimulations of a smaller system (SiO2)0.70(Al2O3)0.13(CaO)0.17 contain-ing 927 atoms are performed here. As discussed in the last section,the system size is not large enough to reproduce the distribution ofcoordination numbers of Ca atoms; however, examples of DFT-MDsimulations of glass materials with this size are still rare.63.3.1. Structural properties of (SiO2 )0.70(Al2O3)0.13(CaO)0.17The initial structure for the present DFT-MD simulation was pre-pared by an NVT-MD simulation with FF-Pedone using the experimen-tal volume. DFT-MD simulations with the fixed volume were performedusing a time step of 2.0 fs and the temperature was controlled at 300 Kusing the stochastic velocity rescaling method [49]. Fig. 6(a) showsthe calculated pressure during the simulation. (See Figs. S1 and S2in Supplementary Material for the temperature and potential energy.)The results show that: (i) the time needed for the equilibration is onlyabout 0.1–0.2 ps, and (ii) the calculated internal pressure in the NVTsimulation with the experimental volume is about −0.7 GPa, closeto 0 GPa. This result suggests that the optimum volume of DFT-MDsimulations with GGA is close to the experimental value. If the celllength is reduced by 0.5%, the calculated stress becomes almost 0GPa, as can be seen by the dotted line in Fig. 6(a). Therefore, DFT-MDsimulations based on GGA can reproduce the experimental volume ordensity, within 1%–2%. The same procedures are followed for vitreousSiO2 and (SiO2)0.5(CaO)0.5 and the calculated stresses are found to benegative and positive, respectively. However, they are both close to 0GPa, about -3 GPa and +0.5 GPa, respectively.Using the trajectory of the DFT-MD simulations of (SiO2)0.70(Al2O3)0.13(CaO)0.17, the structure factor is calculated and shown inFig. 6(b), together with the experimental result of X-ray scattering,measured at Beamline BL04B2 at SPring-8. The calculated result agreeswell with the experimental one in the wide region of wave length 𝑄.The structure factor calculated from a snapshot of NVT-MD simulationswith FF-Pedone also reproduces the experimental result. However, itis found that there are also some important differences between thestructures from the DFT-MD and MD with FF-Pedone.Fig. 7 shows the distribution of interatomic distances and bondangles for the snapshot structures obtained by the MD simulation basedon DFT and FF-Pedone. For the differences in the bond lengths ofSi–O and Al–O, the average lengths by DFT-MD simulations are largerthan those from the MD by FF-Pedone. Note that the optimized latticeparameters of crystalline silica and alumina by GGA with the presentDZP basis set have an error of about 0.5%, as shown in Table 1.However, the difference in the average bond length between FF-MDand DFT-MD is larger than this value, being about 3%. Furthermore,much clearer differences are observed in the distributions of Si–O–Siangle and Ca–O distance. The average Si–O–Si angle is 140◦, in theDFT-MD simulations while it is 160◦, in the MD with FF-Pedone. Thislarge difference in the angle distribution may come from the shorterSi–O distance in the MD simulations using FF-Pedone. It may beworth pointing out that there are no peak structures around 90◦, whichis related to the defect formation of two-membered rings and wasdiscussed in Refs. [50,51].MD simulations of (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 for 𝑋 = Mg, Sr andBa were also performed. Fig. 7(f)–(h) show the distributions of 𝑋–OJournal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 5. (a) Comparison of the atomic forces calculated by DFT (green crosses) and classical force fields (purple open circles) for O atoms in (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg,Ca, Sr, Ba). (b) Difference in the two forces for the bridging oxygen (BO) and non-bridging oxygen (NBO) atoms are plotted using blue open circles and red crosses, respectively.In (b), 200 atoms are selected randomly for both BO and NBO atoms. In both figures, the horizontal axis shows a shorter distance between oxygen and the neighboring siliconatoms.distances obtained by MD simulations using DFT and FF-Pedone. In thecase of 𝑋 = Mg, the distribution is similar for the DFT-MD and MDwith FF-Pedone simulations, while the 𝑋–O distances by the DFT-MDsimulations show wider distributions in other cases (𝑋 = Sr or Ba), ascan be seen in the Ca case.3.3.2. Density of states of (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca,Sr, Ba)The electronic structure of (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 was alsoinvestigated. As shown in Fig. 8, the total density of states (DOS) forthe occupied electrons, ranging from −8 eV to the Fermi level (middleof the band gap), is found to be almost the same for the different 𝑋s.In contrast, the DOS for 𝑋 = Mg at the unoccupied region near theFermi level is completely different from those in other cases. Fig. 8(b)–(e) show the projected DOS for Si, Al, O, and 𝑋 atoms, respectively.The figure shows that the states at the conduction band near the Fermilevel mainly comes from𝑋 for Ca, Sr, and Ba, but no major contributionfrom Mg atoms exist in the same region in the case of 𝑋 = Mg. This7unique electronic property of Mg may be observed in the actual multi-component glass materials; for example, in the absorption spectrum, ifthere are no other peaks from the impurities in this energy region.4. SummaryAlthough molecular dynamics simulations can be a powerful tool toexplain the atomic-scale structures of multi-component glass systems,there are several problems to employ reliable and accurate MD simu-lations of such large and complex materials. The available FFs whichcan treat many elements are limited and, even if such a FF exists,its accuracy is uncertain due to the large variety and complexity ofchemical environments in the multi-component glasses. Although first-principles calculations based on DFT are reliable even for such complexsystems, DFT calculations of multi-component systems are difficultbecause they are large systems containing 1000–5000 atoms or more.Conventional DFT calculations of such large systems are prohibitivelyexpensive.Journal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 6. (a) Calculated pressure during the constant-volume DFT-MD simulations of (SiO2)0.70(Al2O3)0.13(CaO)0.17, with the experimental volume and the reduced volume whosecell parameter is reduced by 0.5%. (b) Structure factor calculated from the snapshots of the DFT-MD simulations, with the experimental data for comparison.Fig. 7. Comparison of distributions in (a) Si–O–Si angle,(b) O–Si–O angle, (c) Si–O distance (d) Al–O distance, (e) Ca–O distance, (f) Mg–O distance, (g)Sr–O distance, and (h)Ba–Odistance by the DFT-MD and FF-MD models. (a), (b), (c), (d) and (e) are calculated from the snapshot structure of (SiO2)0.70(Al2O3)0.13(CaO)0.17, while (f), (g) and (h) are from(SiO2)0.70(Al2O3)0.13(MgO)0.17, (SiO2)0.70(Al2O3)0.13(SrO)0.17, and (SiO2)0.70(Al2O3)0.13(BaO)0.17, respectively.In this study, we used large-scale DFT techniques with the CON-QUEST code to perform a large-scale DFT study of multi-componentglasses. The results demonstrated that accurate DFT calculations ofmulti-component glass systems, (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 =Mg, Ca, Sr, Ba) containing about 5000 atoms are possible using the mul-tisite method, implemented in the CONQUEST code. Using the results ofsuch large-sale DFT calculations, we have evaluated the accuracy of FF-Pedone, FF-Deng, and FF-Wang for vitreous SiO2, (SiO2)0.50(Al2O3)0.50,and (SiO ) (CaO) , and FF-Pedone for (SiO ) (Al O ) (𝑋O)82 0.50 0.50 2 0.70 2 3 0.13 0.17(𝑋 = Mg, Ca, Sr, Ba). By analyzing the trajectories of classical MDsimulations, it was found that the Si–O distance of the non-bridgingoxygen (NBO) atoms is often shorter than 1.55 Å in the simulationsusing FF-Pedone and FF-Deng, while such short Si–O distances werenot observed with FF-Wang. By comparing the forces calculated byFFs and DFT, large deviations from the DFT forces were observed forNBO atoms, specifically for short Si–O bonds, in the case of FF-Pedoneand FF-Deng. However, the accuracy of FF-Wang was found to beindependent of the difference between NBO and BO atoms. DetailedJournal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.Fig. 8. (a) Total density of states (DOS) of (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca, Sr, Ba) and its projected DOS for (b) Si, (c) O, (d) Al, and (e) 𝑋 atoms.).comparison of the three FFs shows that they have clear differenceswhen the calculated forces deviate largely from DFT forces. A similartendency is found for (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 = Mg, Ca, Sr,Ba) with FF-Pedone. The accuracy of the FF is low for NBO atoms,especially when the distance of Si–O is short.DFT-MD simulations of about 1000 atom multi-component glasssystems were also performed and some significant differences in thedistribution of distances and angles between the DFT-MD and FF-MD simulations, specifically around the oxygen atoms. In DFT-MDsimulations, smaller angles (Si–O–Si) and larger bond lengths (Si–O)were observed, and as a result of the cancellation, the distance in theSi atoms is almost the same between DFT-MD and FF-MD simulations,both of which used the experimental volume for the simulation cell.For the 𝑋–O distances, DFT-MD simulations show wider distributionscompared with those FF-MD, except in the case of 𝑋 = Mg, where theagreement between the two MD simulations is much better.It was found that the calculated pressure is very low in the DFT-MD simulations using the experimental volume, and the calculatedstructure factor shows good agreement with the experimental results.Using the snapshot structures obtained by DFT-MD simulations, we alsorevealed the electronic structure of (SiO2)0.70(Al2O3)0.13(𝑋O)0.17 (𝑋 =Mg, Ca, Sr, Ba). The analysis of the total and projected DOS (pDOS)shows that the electronic structure in the occupied region is almost thesame for the different 𝑋s, while the Mg case shows a unique electronicstructure in the unoccupied region near the Fermi level.These results lead to the conclusion that the large-scale DFT meth-ods can be applied to glass systems including multi-component glasses,and can provide various insights, such as the verification of FFs, localstructures around specific elements and electronic structure. Manyanalyses of the mechanical and thermodynamic properties, includingfree energy profiles using the Blue-moon ensemble [52], can be per-formed with large-scale DFT techniques, just as with conventionalDFT calculations. In the near future, large-scale DFT calculations areexpected to make significant contributions to theoretical studies on thedetailed structural and physical properties of multi-component glasses.CRediT authorship contribution statementAtsushi Tanaka: Investigation, Software, Methodology, Formalanalysis, Writing – original draft. Atsuki Saito: Conceptualization,Writing – review & editing, Supervision. Takashi Murata: Supervi-sion, Project administration. Ayako Nakata: Methodology, Software.Tsuyoshi Miyazaki: Conceptualization, Methodology, Writing – origi-nal draft, Supervision.9Declaration of competing interestThe authors declare that they have no known competing finan-cial interests or personal relationships that could have appeared toinfluence the work reported in this paper.Data availabilityData will be made available on request.AcknowledgmentsThis work was supported by World Premier International ResearchCenter Initiative (WPI), MEXT, Japan and JSPS Grant-in-Aid for Trans-formative Research Areas (A) ‘‘Hyper-Ordered Structures Science’’(JP20H05883). Calculations were performed on the Numerical Mate-rials Simulator at NIMS. The synchrotron radiation experiments wereperformed at BL04B2 in SPring-8 with the approval of Japan Syn-chrotron Radiation Research Institute (JASRI) (Proposal No. 2019B2007Appendix A. Supplementary dataSupplementary material related to this article can be found onlineat https://doi.org/10.1016/j.jnoncrysol.2023.122714.References[1] J.E. Shelby, Introduction to glass science and technology, 2005.[2] W.H. Zachariasen, The atomic arrangement in glass, J. Am. Chem. Soc.54 (10) (1932) 3841–3851, http://dx.doi.org/10.1021/ja01349a006, arXiv:10.1021/ja01349a006.[3] S.R. Elliott, Medium-range structural order in covalent amorphous solids, Nature354 (1991) 445–452, http://dx.doi.org/10.1038/354445a0.[4] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, The structural groupsof alkali silicate glasses determined from 29si MAS-nmr, J. Non-Cryst. Solids127 (1) (1991) 53–64, http://dx.doi.org/10.1016/0022-3093(91)90400-Z, URLhttps://www.sciencedirect.com/science/article/pii/002230939190400Z.[5] C.I. Merzbacher, W.B. White, The structure of alkaline earth aluminosilicateglasses as determined by vibrational spectroscopy, J. Non-Cryst. Solids 130(1) (1991) 18–34, http://dx.doi.org/10.1016/0022-3093(91)90152-V, URL https://www.sciencedirect.com/science/article/pii/002230939190152V.[6] A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, A newself-consistent empirical interatomic potential model for oxides, silicates, andsilica-based glasses, J. Phys. Chem. B 110 (2006) 11780–11795.https://doi.org/10.1016/j.jnoncrysol.2023.122714http://refhub.elsevier.com/S0022-3093(23)00579-3/sb1http://dx.doi.org/10.1021/ja01349a006http://arxiv.org/abs/10.1021/ja01349a006http://arxiv.org/abs/10.1021/ja01349a006http://arxiv.org/abs/10.1021/ja01349a006http://dx.doi.org/10.1038/354445a0http://dx.doi.org/10.1016/0022-3093(91)90400-Zhttps://www.sciencedirect.com/science/article/pii/002230939190400Zhttp://dx.doi.org/10.1016/0022-3093(91)90152-Vhttps://www.sciencedirect.com/science/article/pii/002230939190152Vhttps://www.sciencedirect.com/science/article/pii/002230939190152Vhttps://www.sciencedirect.com/science/article/pii/002230939190152Vhttp://refhub.elsevier.com/S0022-3093(23)00579-3/sb6http://refhub.elsevier.com/S0022-3093(23)00579-3/sb6http://refhub.elsevier.com/S0022-3093(23)00579-3/sb6http://refhub.elsevier.com/S0022-3093(23)00579-3/sb6http://refhub.elsevier.com/S0022-3093(23)00579-3/sb6Journal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.[7] L. Deng, J. Du, Development of boron oxide potentials for computersimulations of multicomponent oxide glasses, J. Am. Ceram. Soc. 102(5) (2019) 2482–2505, http://dx.doi.org/10.1111/jace.16082, arXiv:https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.16082, URL https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.16082.[8] M. Wang, N. Anoop Krishnan, B. Wang, M.M. Smedskjaer, J.C. Mauro, M.Bauchy, A new transferable interatomic potential for molecular dynamics sim-ulations of borosilicate glasses, J. Non-Cryst. Solids 498 (2018) 294–304, http://dx.doi.org/10.1016/j.jnoncrysol.2018.04.063, URL https://www.sciencedirect.com/science/article/pii/S0022309318302643.[9] S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, First-principles interatomicpotential of silica applied to molecular dynamics, Phys. Rev. Lett. 61 (1988)869–872, http://dx.doi.org/10.1103/PhysRevLett.61.869, URL https://link.aps.org/doi/10.1103/PhysRevLett.61.869.[10] B.W.H. van Beest, G.J. Kramer, R.A. van Santen, Force fields for silicas andaluminophosphates based on ab initio calculations, Phys. Rev. Lett. 64 (1990)1955–1958, http://dx.doi.org/10.1103/PhysRevLett.64.1955, URL https://link.aps.org/doi/10.1103/PhysRevLett.64.1955.[11] P. Tangney, S. Scandolo, An ab initio parametrized interatomic force field forsilica, J. Chem. Phys. 117 (19) (2002) 8898–8904, http://dx.doi.org/10.1063/1.1513312, arXiv:10.1063/1.1513312.[12] F. Noritake, K. Kawamura, T. Yoshino, E. Takahashi, Molecular dynam-ics simulation and electrical conductivity measurement of na2o*3sio2 meltunder high pressure; relationship between its structure and properties, J.Non-Cryst. Solids 358 (23) (2012) 3109–3118, http://dx.doi.org/10.1016/j.jnoncrysol.2012.08.027, URL https://www.sciencedirect.com/science/article/pii/S0022309312005066.[13] S. Urata, N. Nakamura, T. Tada, H. Hosono, Molecular dynamics study onthe co-doping effect of al2o3 and fluorine to reduce Rayleigh scatteringof silica glass, J. Am. Ceram. Soc. 104 (10) (2021) 5001–5015, http://dx.doi.org/10.1111/jace.17774, arXiv:https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17774, URL https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17774.[14] S. Urata, Modeling short-range and three-membered ring structures in lithiumborosilicate glasses using a machine-learning potential, J. Phys. Chem. C 126(50) (2022/12/22) 21507–21517, http://dx.doi.org/10.1021/acs.jpcc.2c07597.[15] J.C. Fogarty, H.M. Aktulga, A.Y. Grama, A.C.T. van Duin, S.A. Pandit, A reactivemolecular dynamics simulation of the silica-water interface, J. Chem. Phys.132 (17) (2010) 174704, http://dx.doi.org/10.1063/1.3407433, arXiv:10.1063/1.3407433.[16] T. Ohkubo, E. Tsuchida, K. Deguchi, S. Ohki, T. Shimizu, T. Otomo, Y.Iwadate, Insights from ab initio molecular dynamics simulations for a multi-component oxide glass, J. Am. Ceram. Soc. 101 (3) (2018) 1122–1134, http://dx.doi.org/10.1111/jace.15269, arXiv:https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.15269, URL https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.15269.[17] H. Gong, B. Song, Y. Yang, P. Wang, Z. Cao, X. Chen, G. Zhao, S. Peng,Y. Liu, G. Han, Ab initio molecular dynamics simulation of the structuraland electronic properties of aluminoborosilicate glass, J. Am. Ceram. Soc.104 (7) (2021) 3198–3211, http://dx.doi.org/10.1111/jace.17761, arXiv:https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17761, URL https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17761.[18] D.R. Bowler, T. Miyazaki, ( ) Methods in electronic structure calculations,Rep. Progr. Phys. 75 (2012) 36503, http://dx.doi.org/10.1088/0034-4885/75/3/036503.[19] S. Kohara, J. Akola, H. Morita, K. Suzuya, J.K.R. Weber, M.C. Wilding, C.J.Benmore, Relationship between topological order and glass forming ability indensely packed enstatite and forsterite composition glasses, Proc. Natl. Acad. Sci.108 (36) (2011) 14780–14785, http://dx.doi.org/10.1073/pnas.1104692108,arXiv:https://www.pnas.org/doi/pdf/10.1073/pnas.1104692108, URL https://www.pnas.org/doi/abs/10.1073/pnas.1104692108.[20] J. Akola, S. Kohara, K. Ohara, A. Fujiwara, Y. Watanabe, A. Masuno, T. Usuki,T. Kubo, A. Nakahira, K. Nitta, T. Uruga, J.K.R. Weber, C.J. Benmore, Networktopology for the formation of solvated electrons in binary cao-al2o3 compositionglasses, Proc. Natl. Acad. Sci. 110 (25) (2013) 10129–10134, http://dx.doi.org/10.1073/pnas.1300908110, arXiv:https://www.pnas.org/doi/pdf/10.1073/pnas.1300908110, URL https://www.pnas.org/doi/abs/10.1073/pnas.1300908110.[21] K. Konstantinou, P.V. Sushko, D.M. Duffy, Modelling the local atomic structureof molybdenum in nuclear waste glasses with ab initio molecular dynamicssimulations, Phys. Chem. Chem. Phys. 18 (2016) 26125–26132, http://dx.doi.org/10.1039/C6CP03076A.[22] Y. Qian, B. Song, J. Jin, G.I. Prayogo, K. Utimula, K. Nakano, R. Maezono,K. Hongo, G. Zhao, Ab initio molecular dynamics simulation of structuraland elastic properties of SiO2–P2O5–Al2O3–Na2O glass, J. Am. Ceram. Soc.105 (11) (2022) 6604–6615, http://dx.doi.org/10.1111/jace.18614, arXiv:https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.18614, URL https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.18614.10[23] E. Hernández, M.J. Gillan, C.M. Goringe, Linear-scaling density-functional-theorytechnique: The density-matrix approach, Phys. Rev. B 53 (11) (1996) 7147–7157,http://dx.doi.org/10.1103/PhysRevB.53.7147.[24] D.R. Bowler, T. Miyazaki, M.J. Gillan, Recent progress in linear scaling ab initioelectronic structure techniques, J. Phys.: Condens. Matter 14 (2002) 2781–2798.[25] D.R. Bowler, R. Choudhury, M.J. Gillan, T. Miyazaki, Recent progress withlarge-scale ab initio calculations: the conquest code, Phys. Status Solidi (B)243 (5) (2006) 989–1000, http://dx.doi.org/10.1002/pssb.200541386, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200541386, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.200541386.[26] A. Nakata, J.S. Baker, S.Y. Mujahed, J.T.L. Poulton, S. Arapan, J. Lin, Z. Raza,S. Yadav, L. Truflandier, T. Miyazaki, D.R. Bowler, Large scale and linearscaling DFT with the CONQUEST code, J. Chem. Phys. 152 (16) (2020) 164112,http://dx.doi.org/10.1063/5.0005074, arXiv:10.1063/5.0005074.[27] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximationmade simple [Phys. Rev. Lett. 77, 3865 (1996)],, Phys. Rev. Lett. 78 (1997)1396, http://dx.doi.org/10.1103/PhysRevLett.78.1396, URL https://link.aps.org/doi/10.1103/PhysRevLett.78.1396.[28] H. Hay, G. Ferlat, M. Casula, A.P. Seitsonen, F. Mauri, Dispersion effectsin SiO2 polymorphs: An ab initio study, Phys. Rev. B 92 (2015) 144111,http://dx.doi.org/10.1103/PhysRevB.92.144111, URL https://link.aps.org/doi/10.1103/PhysRevB.92.144111.[29] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem. 27 (15) (2006) 1787–1799,http://dx.doi.org/10.1002/jcc.20495.[30] E. Hernández, M.J. Gillan, C.M. Goringe, Basis functions for linear-scaling first-principles calculations, Phys. Rev. B 55 (20) (1997) 13485–13493, http://dx.doi.org/10.1103/PhysRevB.55.13485.[31] D.R. Bowler, T. Miyazaki, A. Nakata, L. Truflandier, URL https://ordern.github.io.[32] D.R. Hamann, URL https://github.com/pipidog/ONCVPSP.[33] D.R. Hamann, Optimized norm-conserving vanderbilt pseudopotentials, Phys.Rev. B 88 (2013) 085117, http://dx.doi.org/10.1103/PhysRevB.88.085117, URLhttps://link.aps.org/doi/10.1103/PhysRevB.88.085117.[34] P. project, URL http://www.pseudo-dojo.org.[35] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X.Gonze, G.-M. Rignanese, The PseudoDojo: Training and grading a 85 elementoptimized norm-conserving pseudopotential table, Comput. Phys. Comm. 226(2018) 39–54, http://dx.doi.org/10.1016/j.cpc.2018.01.012, URL https://www.sciencedirect.com/science/article/pii/S0010465518300250.[36] J.M. Soler, E. Artacho, J.D. Gale, A. García, J. Junquera, P. Ordejón, D. Sánchez-Portal, The SIESTA method for ab initio order-n materials simulation, J. Phys.:Condens. Matter 14 (2002) 2745.[37] A.S. Torralba, M. Todorovic, V. Brázdová, R. Choudhury, T. Miyazaki, M.J.Gillan, D.R. Bowler, Pseudo-atomic orbitals as basis sets for the O(N) DFT codeCONQUEST, J. Phys.: Condens. Matter 20 (29) (2008) 294206, http://dx.doi.org/10.1088/0953-8984/20/29/294206.[38] D.R. Bowler, J.S. Baker, J.T.L. Poulton, S.Y. Mujahed, J. Lin, S. Yadav, Z. Raza,T. Miyazaki, Highly accurate local basis sets for large-scale DFT calculations inconquest, Japan. J. Appl. Phys. 58 (2019) 100503, http://dx.doi.org/10.7567/1347-4065/ab45af.[39] L. Levien, C.T. Prewitt, D.J. Weidner, Structure and elastic properties of quartzat pressure, Am. Mineral. 65 (1980) 920–930.[40] R. Franco, M. Blanco, A. Martín Pendás, E. Francisco, J. Recio, Atomistic simu-lation of the pressure-temperature-volume diagram in 𝛼-Al2O3, Solid State Com-mun. 98 (1) (1996) 41–44, http://dx.doi.org/10.1016/0038-1098(96)00015-4,URL https://www.sciencedirect.com/science/article/pii/0038109896000154.[41] B. Gupta, R. Goyal, Static and thermophysical properties of chalcogenide crystalswith nacl structure, Solid State Commun. 49 (6) (1984) 559–562, http://dx.doi.org/10.1016/0038-1098(84)90191-1, URL https://www.sciencedirect.com/science/article/pii/0038109884901911.[42] A. Nakata, D.R. Bowler, T. Miyazaki, Efficient calculations with multisite localorbitals in a large-scale DFT code CONQUEST, J. Chem. Theory Comput. 10(2014) 4813, http://dx.doi.org/10.1021/ct5004934.[43] A. Nakata, D. Bowler, T. Miyazaki, Optimized multi-site local orbitals in thelarge-scale DFT program CONQUEST, Phys. Chem. Chem. Phys. 17 (2015) 31427,http://dx.doi.org/10.1039/C5CP00934K.[44] M.J. Rayson, P.R. Briddon, Highly efficient method for Kohn-Sham densityfunctional calculations of 500 − −10 000 atom systems, Phys. Rev. B 80 (2009)205104, http://dx.doi.org/10.1103/PhysRevB.80.205104.[45] M. Rayson, Rapid filtration algorithm to construct a minimal basis on the flyfrom a primitive Gaussian basis, Comput. Phys. Comm. 181 (2010) 1051–1056,http://dx.doi.org/10.1016/j.cpc.2010.02.012.[46] G.A. Rosales-Sosa, A. Masuno, Y. Higo, H. Inoue, Crack-resistant Al2O3–SiO2glasses, Sci. Rep. 6 (1) (2016) 1–7, http://dx.doi.org/10.1038/srep23620.http://dx.doi.org/10.1111/jace.16082http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.16082http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.16082http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.16082https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.16082https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.16082https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.16082http://dx.doi.org/10.1016/j.jnoncrysol.2018.04.063http://dx.doi.org/10.1016/j.jnoncrysol.2018.04.063http://dx.doi.org/10.1016/j.jnoncrysol.2018.04.063https://www.sciencedirect.com/science/article/pii/S0022309318302643https://www.sciencedirect.com/science/article/pii/S0022309318302643https://www.sciencedirect.com/science/article/pii/S0022309318302643http://dx.doi.org/10.1103/PhysRevLett.61.869https://link.aps.org/doi/10.1103/PhysRevLett.61.869https://link.aps.org/doi/10.1103/PhysRevLett.61.869https://link.aps.org/doi/10.1103/PhysRevLett.61.869http://dx.doi.org/10.1103/PhysRevLett.64.1955https://link.aps.org/doi/10.1103/PhysRevLett.64.1955https://link.aps.org/doi/10.1103/PhysRevLett.64.1955https://link.aps.org/doi/10.1103/PhysRevLett.64.1955http://dx.doi.org/10.1063/1.1513312http://dx.doi.org/10.1063/1.1513312http://dx.doi.org/10.1063/1.1513312http://arxiv.org/abs/10.1063/1.1513312http://dx.doi.org/10.1016/j.jnoncrysol.2012.08.027http://dx.doi.org/10.1016/j.jnoncrysol.2012.08.027http://dx.doi.org/10.1016/j.jnoncrysol.2012.08.027https://www.sciencedirect.com/science/article/pii/S0022309312005066https://www.sciencedirect.com/science/article/pii/S0022309312005066https://www.sciencedirect.com/science/article/pii/S0022309312005066http://dx.doi.org/10.1111/jace.17774http://dx.doi.org/10.1111/jace.17774http://dx.doi.org/10.1111/jace.17774http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17774http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17774http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17774https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17774https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17774https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17774http://dx.doi.org/10.1021/acs.jpcc.2c07597http://dx.doi.org/10.1063/1.3407433http://arxiv.org/abs/10.1063/1.3407433http://arxiv.org/abs/10.1063/1.3407433http://arxiv.org/abs/10.1063/1.3407433http://dx.doi.org/10.1111/jace.15269http://dx.doi.org/10.1111/jace.15269http://dx.doi.org/10.1111/jace.15269http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.15269http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.15269http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.15269https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.15269https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.15269https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.15269http://dx.doi.org/10.1111/jace.17761http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17761http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17761http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.17761https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17761https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17761https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.17761http://dx.doi.org/10.1088/0034-4885/75/3/036503http://dx.doi.org/10.1088/0034-4885/75/3/036503http://dx.doi.org/10.1088/0034-4885/75/3/036503http://dx.doi.org/10.1073/pnas.1104692108http://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.1104692108https://www.pnas.org/doi/abs/10.1073/pnas.1104692108https://www.pnas.org/doi/abs/10.1073/pnas.1104692108https://www.pnas.org/doi/abs/10.1073/pnas.1104692108http://dx.doi.org/10.1073/pnas.1300908110http://dx.doi.org/10.1073/pnas.1300908110http://dx.doi.org/10.1073/pnas.1300908110http://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.1300908110http://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.1300908110http://arxiv.org/abs/https://www.pnas.org/doi/pdf/10.1073/pnas.1300908110https://www.pnas.org/doi/abs/10.1073/pnas.1300908110http://dx.doi.org/10.1039/C6CP03076Ahttp://dx.doi.org/10.1039/C6CP03076Ahttp://dx.doi.org/10.1039/C6CP03076Ahttp://dx.doi.org/10.1111/jace.18614http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.18614http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.18614http://arxiv.org/abs/https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/jace.18614https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.18614https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.18614https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/jace.18614http://dx.doi.org/10.1103/PhysRevB.53.7147http://refhub.elsevier.com/S0022-3093(23)00579-3/sb24http://refhub.elsevier.com/S0022-3093(23)00579-3/sb24http://refhub.elsevier.com/S0022-3093(23)00579-3/sb24http://dx.doi.org/10.1002/pssb.200541386http://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200541386http://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200541386http://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200541386https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.200541386https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.200541386https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.200541386http://dx.doi.org/10.1063/5.0005074http://arxiv.org/abs/10.1063/5.0005074http://dx.doi.org/10.1103/PhysRevLett.78.1396https://link.aps.org/doi/10.1103/PhysRevLett.78.1396https://link.aps.org/doi/10.1103/PhysRevLett.78.1396https://link.aps.org/doi/10.1103/PhysRevLett.78.1396http://dx.doi.org/10.1103/PhysRevB.92.144111https://link.aps.org/doi/10.1103/PhysRevB.92.144111https://link.aps.org/doi/10.1103/PhysRevB.92.144111https://link.aps.org/doi/10.1103/PhysRevB.92.144111http://dx.doi.org/10.1002/jcc.20495http://dx.doi.org/10.1103/PhysRevB.55.13485http://dx.doi.org/10.1103/PhysRevB.55.13485http://dx.doi.org/10.1103/PhysRevB.55.13485https://ordern.github.iohttps://ordern.github.iohttps://ordern.github.iohttps://github.com/pipidog/ONCVPSPhttp://dx.doi.org/10.1103/PhysRevB.88.085117https://link.aps.org/doi/10.1103/PhysRevB.88.085117http://www.pseudo-dojo.orghttp://dx.doi.org/10.1016/j.cpc.2018.01.012https://www.sciencedirect.com/science/article/pii/S0010465518300250https://www.sciencedirect.com/science/article/pii/S0010465518300250https://www.sciencedirect.com/science/article/pii/S0010465518300250http://refhub.elsevier.com/S0022-3093(23)00579-3/sb36http://refhub.elsevier.com/S0022-3093(23)00579-3/sb36http://refhub.elsevier.com/S0022-3093(23)00579-3/sb36http://refhub.elsevier.com/S0022-3093(23)00579-3/sb36http://refhub.elsevier.com/S0022-3093(23)00579-3/sb36http://dx.doi.org/10.1088/0953-8984/20/29/294206http://dx.doi.org/10.1088/0953-8984/20/29/294206http://dx.doi.org/10.1088/0953-8984/20/29/294206http://dx.doi.org/10.7567/1347-4065/ab45afhttp://dx.doi.org/10.7567/1347-4065/ab45afhttp://dx.doi.org/10.7567/1347-4065/ab45afhttp://refhub.elsevier.com/S0022-3093(23)00579-3/sb39http://refhub.elsevier.com/S0022-3093(23)00579-3/sb39http://refhub.elsevier.com/S0022-3093(23)00579-3/sb39http://dx.doi.org/10.1016/0038-1098(96)00015-4https://www.sciencedirect.com/science/article/pii/0038109896000154http://dx.doi.org/10.1016/0038-1098(84)90191-1http://dx.doi.org/10.1016/0038-1098(84)90191-1http://dx.doi.org/10.1016/0038-1098(84)90191-1https://www.sciencedirect.com/science/article/pii/0038109884901911https://www.sciencedirect.com/science/article/pii/0038109884901911https://www.sciencedirect.com/science/article/pii/0038109884901911http://dx.doi.org/10.1021/ct5004934http://dx.doi.org/10.1039/C5CP00934Khttp://dx.doi.org/10.1103/PhysRevB.80.205104http://dx.doi.org/10.1016/j.cpc.2010.02.012http://dx.doi.org/10.1038/srep23620Journal of Non-Crystalline Solids 625 (2024) 122714A. Tanaka et al.[47] For (SiO2)0.5(Al2O3)0.5, we used 2.78 g/cm3 which was obtained by a linearinterpolation from the experimental densities of the glasses having similarcompositions.[48] V. Dimitrov, T. Komatsu, Electronic polarizability, optical basicity and non-linear optical properties of oxide glasses, J. Non-Cryst. Solids 249 (2)(1999) 160–179, http://dx.doi.org/10.1016/S0022-3093(99)00317-8, URL https://www.sciencedirect.com/science/article/pii/S0022309399003178.[49] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocityrescaling, J. Chem. Phys. 126 (2007) 014101, http://dx.doi.org/10.1063/1.2408420.11[50] M. Kim, K.H. Khoo, J.R. Chelikowsky, Simulating liquid and amorphous silicondioxide using real-space pseudopotentials, Phys. Rev. B 86 (2012) 054104,http://dx.doi.org/10.1103/PhysRevB.86.054104, URL https://link.aps.org/doi/10.1103/PhysRevB.86.054104.[51] K. Shirai, K. Watanabe, H. Momida, S. Hyun, First-principles study on the specificheat jump in the glass transition of silica glass and the prigogine-defay ratio, J.Phys.: Condens. Matter 35 (50) (2023) 505401, http://dx.doi.org/10.1088/1361-648X/acf6ec.[52] T. Hirakawa, D.R. Bowler, T. Miyazaki, Y. Morikawa, L.A. Truflandier, Bluemoon ensemble simulation of aquation free energy profiles applied to monoand bifunctional platinum anticancer drugs, J. Comput. Chem. 41 (22) (2020)1973–1984, http://dx.doi.org/10.1002/jcc.26367.http://dx.doi.org/10.1016/S0022-3093(99)00317-8https://www.sciencedirect.com/science/article/pii/S0022309399003178https://www.sciencedirect.com/science/article/pii/S0022309399003178https://www.sciencedirect.com/science/article/pii/S0022309399003178http://dx.doi.org/10.1063/1.2408420http://dx.doi.org/10.1063/1.2408420http://dx.doi.org/10.1063/1.2408420http://dx.doi.org/10.1103/PhysRevB.86.054104https://link.aps.org/doi/10.1103/PhysRevB.86.054104https://link.aps.org/doi/10.1103/PhysRevB.86.054104https://link.aps.org/doi/10.1103/PhysRevB.86.054104http://dx.doi.org/10.1088/1361-648X/acf6echttp://dx.doi.org/10.1088/1361-648X/acf6echttp://dx.doi.org/10.1088/1361-648X/acf6echttp://dx.doi.org/10.1002/jcc.26367 Large-scale DFT calculations of multi-component glass systems (SiO2)0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) : Accuracy of classical force fields Introduction Computational methods Classical Force Fields Large-scale DFT methods Local basis sets and support functions Multisite method for large-scale DFT calculations Modeling of the glass structures by melt-quench process Results Benchmark tests for classical force fields in vitreous silica and binary glasses Large scale DFT calculations of the multicomponent glass systems, (SiO2)0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) DFT-MD simulations of (SiO2 )0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) Structural properties of (SiO2 )0.70(Al2O3)0.13(CaO)0.17 Density of States of (SiO2)0.70(Al2O3)0.13(XO)0.17 (X = Mg, Ca, Sr, Ba) Summary CRediT authorship contribution statement Declaration of competing interest Data availability Acknowledgments Appendix A. Supplementary data References