Special Glass Structures for First Principles Studies of Bulk Metallic Glasses
Siya Zhu, Jan Schroers, Stefano Curtarolo, Hagen Eckert, Axel van de, Walle

TL;DR
This paper introduces a method to design small-cell structures that accurately replicate the local geometric descriptors of bulk metallic glasses, enabling efficient property calculations with ab initio methods.
Contribution
It presents a novel approach combining MD and RMC to optimize small-cell models for bulk metallic glasses, reducing computational costs.
Findings
Successfully reproduces local geometric descriptors of large-cell simulations
Enables accurate property calculations with smaller, computationally manageable models
Implemented in the MAST software package for broader use
Abstract
The atomic-level structure of bulk metallic glasses is a key determinant of their properties. An accurate representation of amorphous systems in computational studies has traditionally required large supercells that are unfortunately computationally demanding to handle using the most accurate ab initio calculations. To address this, we propose to specifically design small-cell structures that best reproduce the local geometric descriptors (e.g., pairwise distances or bond angle distributions) of a large-cell simulation. We rely on molecular dynamics (MD) driven by empirical potentials to generate the target descriptors, while we use reverse Monte Carlo (RMC) methods to optimize the small-cell structure. The latter can then be used to determine mechanical and electronic properties using more accurate electronic structure calculations. The method is implemented in the Metallic Amorphous…
| C11(GPa) | C12(GPa) | K(GPa) | G(GPa) | E (GPa) | ||
|---|---|---|---|---|---|---|
| 1 | 141.34 | 89.77 | 106.96 | 25.78 | 71.60 | 0.388 |
| 2 | 142.11 | 89.23 | 106.86 | 26.44 | 73.27 | 0.386 |
| 3 | 142.30 | 92.63 | 109.19 | 24.83 | 69.25 | 0.394 |
| 4 | 142.13 | 92.88 | 109.30 | 24.63 | 68.72 | 0.395 |
| K(GPa) | G(GPa) | E (GPa) | ||||
|---|---|---|---|---|---|---|
| Experiment [37] | 101.2 | 31.3 | 84.0 | 0.35 | ||
| MAST(32 atoms, average) | 108.08 | 25.42 | 70.7 | 0.391 | 0.998 | 0.998 |
| MAST(12 atoms) | 111.9 | 24.5 | 68.6 | 0.398 | 0.987 | 0.986 |
| EAM-MD (16,000 atoms) | 112.3 | 38.9 | 104.5 | 0.344 | 0.999 | 0.998 |
| EAM-MD + ab initio (32 atoms) | 107.3 | 15.9 | 45.3 | 0.430 | 0.991 | 0.990 |
| EAM-MD (32 atoms) | 110.6 | 40.7 | 108.7 | 0.336 | 0.926 | 0.911 |
| K(GPa,cal) | K(GPa,exp) | G(GPa,cal) | G(GPa,exp) | E (GPa,cal) | E(GPa,exp) | (cal) | (exp) | |
|---|---|---|---|---|---|---|---|---|
| Zr36Cu64 [28, 38] | 104.3 | 34 | 92.3 | 0.352 | ||||
| Zr43.7Cu56.3 | 112.9 | 25.3 | 70.6 | 0.396 | ||||
| Zr46.9Cu53.1 | 112.8 | 20.2 | 57.1 | 0.416 | ||||
| Zr50Cu [37] | 108.8 | 101.2 | 25.4 | 31.3 | 70.7 | 84.0 | 0.391 | 0.35 |
| Zr53.1Cu46.9 | 106.9 | 18.5 | 52.6 | 0.418 | ||||
| Zr56.3Cu43.7 | 107.8 | 21.1 | 59.5 | 0.408 | ||||
| Zr54Cu46 [39, 38] | 128.5 | 30.0 | 83.5 | 0.391 | ||||
| Zr47.5Cu47.5Al5 [37] | 113.7 | 33.0 | 90.1 | 0.365 | ||||
| Zr47.5Cu47.5Al5 [40] | - | - | 60† | - | ||||
| Zr47.5Cu47.5Al5 [41] | 1171 | 321 | 881 | - | ||||
| Zr47.5Cu47.5Al5 [42] | - | - | 875 | - | ||||
| Zr47Cu47Al6 [42] | - | - | 845 | - | ||||
| Zr46.9Cu46.9Al6.2 | 108.2 | 22.3 | 62.6 | 0.40 | ||||
| Zr45Cu45Al10 [43] | - | - | 60† | - | ||||
| Zr45Cu45Al10 [40] | - | - | 80† | - |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsTheoretical and Computational Physics · Metallic Glasses and Amorphous Alloys · nanoparticles nucleation surface interactions
Special Glass Structures for First Principles Studies
of Bulk Metallic Glasses
Siya Zhu
School of Engineering, Brown University, Providence, RI, 02912, USA
Jan Schroers
Department of Mechanical Engineering and Materials Science, Yale University, New Haven, CT, 06511, USA
Stefano Curtarolo
Department of Mechanical Engineering and Materials Science and Center for Autonomous Materials Design, Duke University, Durham, NC, 27708, USA
Hagen Eckert
Department of Mechanical Engineering and Materials Science and Center for Autonomous Materials Design, Duke University, Durham, NC, 27708, USA
Axel van de Walle
School of Engineering, Brown University, Providence, RI 02912, USA
Abstract
The atomic-level structure of bulk metallic glasses is a key determinant of their properties. An accurate representation of amorphous systems in computational studies has traditionally required large supercells that are unfortunately computationally demanding to handle using the most accurate ab initio calculations. To address this, we propose to specifically design small-cell structures that best reproduce the local geometric descriptors (e.g., pairwise distances or bond angle distributions) of a large-cell simulation. We rely on molecular dynamics (MD) driven by empirical potentials to generate the target descriptors, while we use reverse Monte Carlo (RMC) methods to optimize the small-cell structure. The latter can then be used to determine mechanical and electronic properties using more accurate electronic structure calculations. The method is implemented in the Metallic Amorphous Structures Toolkit (MAST) software package.
††preprint: APS/123-QED
I Introduction
For decades, bulk metallic glasses (BMGs) have been investigated for their excellent strength and fracture toughness [1, 2, 3, 4, 5, 6, 7]. Traditional amorphous alloys require high quenching rates to avoid crystallization, which restricts the size of powders and films [8]. In contrast, the critical cooling rates of BMGs are much lower, allowing the size to be greater than in the smallest dimension, and promising the application as structural materials. [9, 10, 11]
The atomic structures of BMGs are of prime importance for predicting and explaining their macroscopic properties [12, 13, 14]. Experimental characterizations with X-ray diffraction (XRD), small-angle neutron scattering(SANS), and transmission electron microscopy (TEM) help us understand some structural properties [15], yet they are insufficient to fully reconstruct the atomic packing of BMGs. Molecular dynamics (MD) is often used in the investigation of atomic structures of BMGs [16]. Traditional embedded atom method (EAM) MD can efficiently simulate the cooling process over a time scale of for systems of over 10,000 atoms [13]. While this is typically sufficient to reproduce the general atomic structure of BMG, EAM-MD exhibits two main shortcomings: First, due to the limited accuracy of EAM, the finer details of the structure, such as specific bond lengths or bond angles, may not be accurately reproduced. As a result, the accuracy of some physical properties, particularly viscosity and plastic mechanical properties, as well as normal mode frequencies, or even material density, may suffer. Second, the structures obtained from EAM-MD cannot, due to their sizes, be directly used for further, more accurate analysis using first-principles calculations. This, therefore, precludes the calculations of electronic properties as well as more accurate calculations of mechanical properties.
To address this, it has been proposed [12] to use first-principles methods for both the simulation of the quenching process and the calculation of the properties of interest. While first-principles methods guarantee a higher accuracy in calculating the energies and forces, the system size is typically limited to atoms and the time scale to picoseconds. As a result, the large statistical fluctuation and periodic boundary conditions in such a small system size as well as the high quenching rate due to the short simulation time can significantly affect the structure of the simulated system.
In an effort to leverage the advantages of both EAM-MD and fully first-principles methods, we propose to combine these methods with the help of Reverse Monte Carlo (RMC). The idea of RMC has often been used to infer the atomic structures of BMG and other disordered systems [17, 18, 19] from indirect experimental characterization techniques. For instance, X-ray diffraction, neutron diffraction, extended X-ray absorption fine structure (EXAFS), etc., can provide structure factors, radial distribution functions, or atomic coordinations. RMC modeling, a variation of the standard Metropolis Monte Carlo method, can then be used to generate structural model (a cell or an ensemble of atoms) by minimizing an error function (also called loss function) quantifying the consistency of the generated structure with the experimental data, while subjecting to a set of constraints like concentrations or cell shape. Various programs implementing this algorithm, such as fullrmc [20], RMCProfile [21], and RMC++ [22] are widely used to study the atomic modeling of liquids, glasses, polymers and other amorphous materials.
Here we use RMC in differently. Instead of using experimental data as input, we perform RMC based on structures from EAM-MD results in order to build representative atomic structures that are small enough for DFT calculations, while maintaining the statistical geometrical properties of the large structure from EAM-MD.
We have developed a C++ software package — Metallic Amorphous Structures Toolkit (MAST), as a part of Alloy Theoretic Automated Toolkit (ATAT) [23, 24]. MAST includes RMC codes for generating random structures suitable for MD and MC run, as well as scripts to automate the entire process, including MD, MC and ab initio calculations. It also contains scripts for calculating mechanical properties such as elastic constants.
We first generate a typical BMG atomic configuration (with over 10,000 atoms) using EAM-MD. Then, we use the RMC algorithm to construct a small periodic system of 32 atoms consistent with the MD result in terms of statistical properties such as pair distribution functions and bond angles. We give the name Special Glass Structure (SGS) to such atomic structure, paraphrasing the concept of Special Quasirandom Structures (SQS) [25, 26, 27] for crystalline structures. Once an optimal small-cell representation has been determined, ab initio calculations for structural optimization and mechanical properties computations are performed.
In this work, we study, as an example, Zr-Cu and Zr-Cu-Al based BMGs over a range of compositions, since these systems represent the transition metal BMG former and are the basis of various important BMGs with attractive properties [28, 29]. Atomic structures and mechanical properties are calculated and compared with experimental data and other theoretical predictions.
II Methods
II.1 Error Functions
To measure the structural differences between a large target structure with over 10,000 atoms and an SGS with 32 atoms, we devise an error function containing various kinds of commonly used structural information.
This information includes the pair distribution functions (also known as the partial radial distribution functions, PDF), defined as
[TABLE]
where is the atomic number density, is the average number of atoms type within a distance of from an atom of type , and is the concentration of atom type in the structure. To make it symmetric to and , it can also be written as
[TABLE]
where is the volume, is the total number of atoms, is the total number of and atoms pairs within distance to , and are the total numbers of atoms type and , respectively. A traditional way to compare the PDFs is to obtain the structure factors with a Fourier transformation on the PDFs, sum over all partial structure factors for the total structure factors, and perform a test with the one from experiments like XRD. Since we only compare between both atomic structures instead of experimental data, we can compute the difference between two PDFs directly. In order to get a smooth PDF without sharp peaks for small structures or noises for large systems, we perform kernel smoothing [30]:
[TABLE]
where is a symmetric kernel function integrating to one and is a user-specified smoothing parameter.
Then, the error function for the PDFs is defined as
[TABLE]
where is a cutoff distance (usually set as the distance of the 2nd or 3rd nearest neighbors), is the PDF of a target structure. We include a factor to give higher weights to shorter distances, as the chemical environment closer to an atom is more important than farther away. Notice that, since the PDF is symmetric in and , we do not double count the term in the sum. For example, for a Zr-Cu system, we only sum over three terms of partial PDFs: Zr-Zr, Zr-Cu, and Cu-Cu, but not Cu-Zr. For a ternary system like Zr-Cu-Al, we sum over 6 terms in .
For bond angles, we define a statistical quantity for the angles between two bonds connecting an atom to two nearest neighbors site occupied by species and :
[TABLE]
where is the number of angles forming with bond - and -, within the range of to . The bond length should be less than a cutoff length, usually set between the distance of the 1st and 2nd nearest neighbours. We also kernel smooth the function . The error function for angles is defined as
[TABLE]
Notice that and are equivalent and is not, therefore, there are 6 terms in the sum for a binary system and 18 terms for a ternary system.
For a total error function, we weigh the errors from each aspect:
[TABLE]
A basic principle for choosing is to make every term of the same order of magnitude, so that they are all taken into account in the Monte Carlo simulations.
II.2 RMC process
With an initial guess, which can be either a MD simulated structure, a previous MC result, or a random structure, we do an RMC process:
- i.
Randomly pick an operation from: (1) move an atom by a random amount; (2) swap positions of two atoms; (3) reshape the supercell by a random amount. The probabilities of each operation are chosen for better convergence. A suggested improvement is to set the initial structure as a cubic cell with volume based on the density of target structure, so that we do not need to reshape the cell (which is a default option in our code). Reshaping cells is sometimes detrimental to the convergence of the method and is not recommended;
- ii.
Calculate the error function for the new structure with the target structure;
- iii.
Compare with the error function of the old structure. If , we always accept the operation. Otherwise, we accept the operation with a probability , where is a “temperature” factor greater than 0. The larger is, the more likely to get over barriers and find the global minimum, while a smaller speeds up the convergence to a local minimum.
- iv.
Go back to the first step, until we reach the we set, or the maximum number of Monte Carlo steps.
II.3 Generation of atomic modeling
As a specific example, to generate a 32-atom SGS with concentration Zr50Cu50, we follow this process:
- i.
Use MAST to generate a random structure with 8000 Zr atoms and 8000 Cu atoms;
- ii.
Use LAMMPS [31, 32] to run an EAM-MD simulation. The EAM potentials in this work are from the group of Sheng [13]. We first do a NVT and a NPT run at , and quenching to at K/s with NPT ensemble. Then a NPT at is performed for a stable structure;
- iii.
Use MAST for an RMC process from a cell with 16 Zr atoms and 16 Cu atoms randomly arranged to generate an SGS;
- iv.
Use LAMMPS to do a structural optimization for the SGS from step(iii) with command “minimize”, then a NPT MD at . This is because RMC may sometimes give unreasonable structures (with two atoms too close to each other). An optimization and MD run with EAM help us avoid this situation;
- v.
Another RMC run on the result of step (iv). More loops of MD and RMC (step (iv) and step (v))can be performed for a better convergence. In this work, we do not perform extra loops of MD and RMC since it already gives very good fitting.
- vi.
Use Vienna ab initio Simulation Package (VASP) [33, 34, 35] to do an ab initio optimization and get the final structure for Zr16Cu16. We use the Perdew-Burke-Ernzerhof (PBE) exchange and correlation functional at level of the generalized gradient approximation(GGA) [36] . A cutoff energy of 600eV is set to eliminate the effect of Pulay stress in structural optimization with ISIF=3. A K-points mesh is set.
II.4 Calculation of Elastic Moduli
As an isotropic material, a BMG only has two independent elastic constants: and . To calculate them, we apply strains on the relaxed structure, optimize with VASP (ISIF = 2), calculate the energy and fit the strain-energy relationship. A planewave energy cutoff of 300eV is used in the mechanical properties calculations and other settings are the same as previously stated. We apply uniaxial strain from to on , and directions (the strain-energy curves on each direction are almost the same, which will be shown in the next section, as an evidence of isotropy of our structure). can be derived from
[TABLE]
Then we add a bulk strain from to on each direction and plot the strain-energy curve. Then can be calculated as
[TABLE]
The bulk modulus is
[TABLE]
and shear modulus is
[TABLE]
III Results
In this section, we first use an example of Zr50Cu50 BMG to illustrate our results. First, we perform a MD simulation to generate a structure with 8000 Zr and 8000 Cu atoms, as shown in Figure 1a. Pair distribution functions and angle distribution functions are calculated and shown in Figure 1b and 1c.
Following the steps in Section II.3, we generate the atomic structures with steps (iii) to (v), calculate pair and angle distribution functions as shown in Figure 2. The RMC steps sometimes generate structures with nonzero PDF at very short distances, thus implying that some atoms are nonphysically close to each other. This can be easily remedied either by performing relaxations and MD of the atoms using the EAM model, as illustrated in Figure 2d-2f. Another RMC process reducing the error function gives an SGS Figure 2g with good fitting, as shown in Figure 2h and 2i. The SGS is now physically reasonable and successfully represents the structural information of the structure in Figure 1a.
After the structural optimization with VASP, we apply uniaxial strain and triaxial strain to our SGS, optimize the structure under strain, and calculate the energy, as shown in Figure 3. Two things need to be noted here: (i) we add the strain from to , which might exceed the range of elastic deformation of BMG. However, in this case, the SGS is too small to perform any plastic deformation or fracture, and the energies still fit well with large strain. Therefore we keep the data of large strain up to 0.03 in this example. The strain range can be set manually in MAST; (ii) all ab initio calculations of energy here are at the temperature of . That would bring about some error in the elastic modulus compared to experimental data at room temperature. Nevertheless, the differences between mechanical properties at and are insignificant. Figure 3b shows the energy under uniaxial strain applied to , and directions, respectively. The strain-energy curves of three directions almost overlap, manifesting that our SGS is “isotropic” in some sense, which provides additional evidence of its validity.
In Table 1, we perform a robustness test of our method on 32 atoms structure of Zr-Cu BMG at concentration . With four different initial structures in step (iii) in Section II.3, we get four different SGSs and calculate their mechanical properties. The results show our Monte Carlo process converges and is robust and reliable.
In Table 2, we compare the mechanical properties of Zr50Cu50 from experiments, our method, and some other methods to show the advantage of our method. For the SGS in Figure 2(c) with 32 atoms, the differences in bulk modulus, shear modulus, Young’s modulus and Poisson’s ratios are less than 15% (labeled MAST (32 atoms) in Table 2). For comparison, we apply strain on the structure with 16,000 atoms in Figure 1a, run MD with EAM potentials for energies, and calculate mechanical properties(labeled as EAM-MD (16,000 atoms) in Table 2). The results from our process MAST agree with experimental data a little bit better — both methods overestimate the bulk modulus K, while MAST underestimates shear and Young’s modulus and EAM-MD overestimates them, compared with experiments. To show the beneficial contribution of our Monte Carlo process in MAST, we also show computations for the glass structure with 32 atoms, generated with EAM-MD only. We calculate the mechanical properties with EAM-MD itself(labeled as EAM-MD (32 atoms) in Table 2) and ab initio (labeled as EAM-MD + ab initio (32 atoms) in Table 2). Both of them show poor agreement with experiments. An SGS with only 12 atoms (Zr6Cu6) is also generated with MAST(labeled as MAST(12 atoms) in Table 2). The mechanical properties of it are close to MAST 32 atoms SGS, and are even better than EAM-MD structures with 32 atoms, which shows the robustness and accuracy of our method.
Similarly, our method can be applied to BMG systems with different concentrations, or even more than two elements. We show some calculation results for ZrCu and ZrCuAl BMG in Table 3, together with some experimental data as reference. For most of the concentrations, our results agree well with experiments on the bulk moduli, but the levels of agreement decrease for shear moduli and Young’s moduli. This originates, in large part, from the fact that ab initio bulk modulus calculations exhibit a higher signal-to-noise ratio, because a triaxial strain induces a larger change in energy than a uniaxial strain and only the bulk modulus depends solely on triaxial strain data (c.f. Eq. 10). As seen in Figure 3, the overall energy changes associated with triaxial strains are about 10 times bigger than those associated with uniaxial strains. In spite of this, the SGS approach still provides sufficient accuracy to investigate the compositional-dependence of elastic properties and, in this case, provides corroboration of the relatively small sensitivity to composition found in experiments.
The usefulness of SGS is not limited to mechanical properties. For instance, they make it possible to analyze detailed electronic structure of BMGs, which would be impossible with EAM alone and intractable with large supercells using ab initio methods. In Figure 4a we show the calculated DOS for different compositions of Zr-Cu glasses. This type of analysis enables us, for instance, to assess whether this BMG is well described by a rigid-band approximation. In the case of Zr-Cu, we can actually determine that shifts in the DOS as a function of composition, relative to the Fermi level, are not uniform over all energies. For instance, around 3eV, there is a larger shift between Zr14Cu18 and Zr15Cu17 composition than between Zr15Cu17 and Zr18Cu14 composition, while to opposite is true around eV.
These effects would be difficult to isolate using other methods. For instance, in Figure 4b, we make a comparison between the DOS obtained via an SGS from MAST and the one obtained via EAM-MD only, showing that the electronic properties would likely be incorrect without the Monte Carlo optimization process. The two DOS differ by shifts of up to 0.2 eV on the energy scale, which is large enough to mask the effects identified in Figure 4a.
IV Conclusion
In this work, we develop a method to generate small-cell approximations to amorphous structures suitable for efficient ab initio calculations. These so-called Special Glass Structures (SGS) are constructed by matching geometric descriptors from large-cell empirical potential MD simulations using reverse Monte Carlo. The method is benchmarked by comparing the predicted mechanical properties of common Bulk Metallic Glasses to corresponding experimental results.
More generally, our method can be used for the study of any amorphous system, and is applicable to the calculation of properties besides the elastic constants, such as thermodynamics properties like enthalpy and vibrational entropy with a higher accuracy than classical potential-based MD. Larger SGS ( atoms) should also apply to the study of kinetic processes, such as viscous flow or plastic deformation processes via shear band evolution in BMGs [44]. While, we leave the demonstration of these capabilities for subsequent work, it should be clear that the SGS method opens the way to the systematic study of BMG properties via high-throughput ab initio methods over a broad range of chemistries [24, 45].
Acknowledgements.
This work is supported by Office of Naval Research grant N00014-20-1-2225. Computational resources were provided by (i) the Center for Computation and Visualization at Brown University, (ii) the Extreme Science and Engineering Discovery Environment (XSEDE) through allocation TG-DMR050013N, which is supported by National Science Foundation Grant No. ACI-1548562 and (iii) the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program through allocation DMR010001, which is supported by National Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Christopher A Schuh, Todd C Hufnagel, and Upadrasta Ramamurty. Mechanical behavior of amorphous alloys. Acta Materialia , 55(12):4067–4109, 2007.
- 2[2] A Reza Yavari, JJ Lewandowski, and J Eckert. Mechanical properties of bulk metallic glasses. Mrs Bulletin , 32(8):635–638, 2007.
- 3[3] J Eckert, J Das, S Pauly, and C Duhamel. Mechanical properties of bulk metallic glasses and composites. Journal of materials research , 22(2):285–301, 2007.
- 4[4] A Lindsay Greer and Evan Ma. Bulk metallic glasses: at the cutting edge of metals research. MRS bulletin , 32(8):611–619, 2007.
- 5[5] MF Ashby and A Lindsay Greer. Metallic glasses as structural materials. Scripta Materialia , 54(3):321–326, 2006.
- 6[6] Ling Shao, Jittisa Ketkaew, Pan Gong, Shaofan Zhao, Sungwoo Sohn, Punnathat Bordeenithikasem, Amit Datye, Rodrigo Miguel Ojeda Mota, Naijia Liu, Sebastian Alexander Kube, et al. Effect of chemical composition on the fracture toughness of bulk metallic glasses. Materialia , 12:100828, 2020.
- 7[7] Henry J Neilson, Alex S Petersen, Andrew M Cheung, S Joseph Poon, Gary J Shiflet, Mike Widom, and John J Lewandowski. Weibull modulus of hardness, bend strength, and tensile strength of ni- ta- co- x metallic glass ribbons. Materials Science and Engineering: A , 634:176–182, 2015.
- 8[8] Jan Schroers. Processing of bulk metallic glass. Advanced materials , 22(14):1566–1597, 2010.
