Elinvar effect in $\beta-$Ti simulated by on-the-fly trained moment tensor potential

A combination of quantum mechanics calculations with machine learning (ML) techniques can lead to a paradigm shift in our ability to predict materials properties from first principles. Here we show that on-the-fly training of an interatomic potential described through moment tensors provides the same accuracy as state-of-the-art {\it ab inito} molecular dynamics in predicting high-temperature elastic properties of materials with two orders of magnitude less computational effort. Using the technique, we investigate high-temperature bcc phase of titanium and predict very weak, Elinvar, temperature dependence of its elastic moduli, similar to the behavior of the so-called GUM Ti-based alloys [T. Sato {\ it et al.}, Science {\bf 300}, 464 (2003)]. Given the fact that GUM alloys have complex chemical compositions and operate at room temperature, Elinvar properties of elemental bcc-Ti observed in the wide temperature interval 1100--1700 K is unique.


2
Theoretical predictions of materials properties have become increasingly more reliable with the development of ab inito molecular dynamics (AIMD), which brings simulations conditions in accord with the conditions at which materials exist in nature or operate in devices. Power of the approach has been demonstrated in numerous applications, ranging from studies of materials at the Earth's core conditions [1] to investigations of alloys for hardcoatings applications [2]. Exploring finite-temperature effects via MD simulations allows one to uncover exceptional effects at finite temperature, such as a collective superionic-like diffusion of atoms in Fe [3] and Ti [4]. However, AIMD calculations are too time-consuming, and this greatly limits their use in practice, despite increasing demands in providing reliable theoretical data for the broad range of applications in physics, chemistry, geosciences, materials science, and biology.
In particular, first-principles simulations of elastic properties of materials at finite temperature are becoming broadly relevant for fundamental research [5] and knowledge-based materials design [6,7]. However, their accurate predictions represent an enormous computational task [8,9]. One has to perform multiple molecular dynamics calculations for relatively large, with hundreds of atoms, simulation cells applying different distortions and explore the stress-strain relations as a function of temperature. Moreover, going from a description of thermodynamic properties of materials at finite temperature to simulations of their elastic constants the accuracy of calculations has to be improved by an order of magnitude, from about 10 meV for potential energy differences to 1 meV. This puts tough requirements on the technical parameters used in AIMD runs, as has been demonstrated, e.g. in studies of the materials with strong elastic softening and low thermal expansion [10].
A traditional way of increasing efficiency of molecular dynamics simulations is to employ classical interatomic potentials fitted to experiment and/or first-principles calculations. The 3 drawback of utilizing the (semi-)empirical interatomic potentials (or force-fields) with MD simulations instead of performing AIMD is that off-the-shelf potentials typically do not have sufficient transferability, and therefore their use for a descriptions of atomic configurations outside of the phase-space used for the original fitting may be problematic. To overcome the issue, the so-called learn-on-the-fly strategy has been suggested [11]. For example, augmentation of the classical potentials via on-the-fly quantum mechanical calculations has been utilized to model brittle failure in silicon [12]. However, the flexibility of such classical potentials is still limited by their functional form.
Machine-learning (ML) of interatomic potentials (MLIP) have been put forward recently as an alternative and highly promising approach what should combine high efficiency and sufficient accuracy in solving relevant physical problems [13][14][15][16][17]. Unlike empirical potentials with a fixed functional form, their ML counterparts aim at approximating an arbitrary interatomic interaction, e.g., by expanding it in a complete set of basis functions whose arguments are local descriptors of atomic environments. In particular, moment tensors have been introduced as descriptors [18]. The corresponding ML many-body interatomic potentials have been called MTPs. They have been successfully applied in zero temperature crystal structure prediction [19], predicting thermal conductivity of compounds [20] or computing anharmonic free energy of refractory high entropy alloys [21]. In the applications, the MTPs are fitted on-the-fly to high-accuracy density functional theory (DFT) calculations and provide ab initio level of accuracy in the description of thermodynamic properties of materials using several orders of magnitude less computational resources than conventional AIMD simulations. However, as has been specified above, calculations of elastic properties of materials typically require order of magnitude more accurate calculations. To the best of our knowledge, applications of MLIP have not been explored for such challenging tasks yet. In this study we employ MTPs in simulations of the elastic constants C ij of bcc-Ti (β-Ti) in the broad temperature interval, from 900 K to 1700 K. Ti is the base elements for a broad family of technologically relevant alloys for airspace [22,23], nuclear energy [24], hard-coatings [25] and biomedical [26,27] applications. From the fundamental point of view, bcc, or β−phase of Ti is highly challenging for theoretical simulations because of its dynamical instability at zero temperature [28]. This means that bcc-Ti cannot exist at low temperature, but it becomes stable not only dynamically, but also thermodynamically at temperatures above 1155 K (at ambient pressure). Therefore, strictly speaking, any theoretical consideration of bcc Ti and Ti-rich alloys requires finite-temperature simulations [29].
We demonstrate that our technique based on the use of MTPs predicts high-temperature elastic properties of bcc-Ti with the same accuracy as state-of-the-art AIMD simulations, though with two orders of magnitude less computational efforts. Importantly, in our simulations we observe very negligible temperature dependence of elastic moduli of β−Ti in the wide temperature interval 1100-1700 K. The Elinvar effect in bcc-Ti uncovered in our simulations is similar to the elastic behavior of the so-called GUM alloy [10,[30][31][32] with an important difference: while the latter have complex chemical compositions and operate at room temperature, the former is a pure elemental metal showing remarkable properties at very high temperatures.
In this work, β-Ti has been modeled with a supercell of 128 atoms. The exchange correlation potential in all our VASP [33] calculations has been approximated using PBE functional [34]. For AIMD simulations we applied the energy cutoff 460 eV and a Γ-point centered (2×2×2) k-point mesh. The elastic constants between 900 and 1700 K have been calculated from the stress-strain relations applying lattice distortions ±1% and ±2% and 5 the NVT corresponding stress-fluctuation formula. We performed 10 000 AIMD steps with each distorted supercells at each temperatures using the time step of 1 fs, see Ref. [35]. The averaging for MTP was done on 100 Metropolis-adjusted Langevin algorithm (MALA) [36] trajectories with 10 5 steps for each distortions.
In our MLIP approach we have used MTP built from polynomial-like basis functions of degree 16 or less, see Ref. [19], which resulted in about 120 parameters to fit. Our twostage active-learning process is illustrated in Fig. 1. The initial training, see Fig. 1 a), set contained six (three temperatures and two different strains) ideal atomic configurations with the experimental bcc lattice parameters [37][38][39] at the corresponding temperature. For the training set we sampled configurations at 900, 1300 and 1700 K and applied ± 2% strains on the simulation box. We performed the corresponding six MD simulations in parallel using LAMMPS [40]. The extrapolation grade of the MTP in each MD at each time step was calculated by the MaxVol algorithm [41]. The arithmetic average of the (six) extrapolation rates is shown in Fig. 1 An MD was stopped at a configuration whose the grade exceeded the extrapolation threshold, Fig. 1 a), and static DFT energy calculation was performed on that extrapolative configuration. Then the training set was expanded the extrapolative configurations (we used at most six configurations) and the MTP was retrained. This process was repeated until the extrapolation rate (inverse of the reliable simulation time) reached zero for each MD.  Fig. 1 b). Thus, we estimate the computational efforts required for the full determination of our MTP potential to correspond to ∼200 accurate DFT calculations.
Let us next verify the accuracy of our potentials in a comparison with available experiemtal information. with just a few hundred static DFT calculations they offer excellent accuracy despite the 7 small expansion coefficient and dynamical instability of bcc-Ti at zero temperature, both of which have adverse effects on the statistical error of the AIMD method.
Next, we utilized our MTP to obtain phonon dispersion relations in β-Ti calculating velocity autocorrelation function and using normal-mode-decomposition implemented into DynaPhoPy [43]. Figure 3 shows the derived phonon dispersions and linewidths extracted from the force constants of a (6 × 6 × 6) supercell and velocity autocorrelation functions calculated with MD simulations of a a (12 × 12 × 12) supercell at 900 and 1500 K.
Quasiharmonic phonons are calculated using the finite displacement approach. Assumption of a harmonic Hamiltonian with force constants derived from the obtained MTP combined with the thermal expansion of the lattice shown in Fig. 2 results in the well-known strong dynamical instability of bcc Ti [44] around N-point and between P an H points.
The corresponding large phonon linewidths (or short lifetimes) are shown with the shaded areas in Fig. 3. It is known that the N-point instability promotes the bcc-hcp transition [45] while the other one is connected to the formation of the ω phase. On the contrary, the phonon dispersions relations calculated from full MD with MTP potentials show the dynamical stability of bcc Ti at high temperature, and are in agreement with experiment and with calculated data from literature [4,44]. Comparison of the results calculated at 900 and 1500 K points to an interesting effect: the broadening decreases with temperature, which is quite unusual. The effect might be related to the fact that the bcc phase of Ti has barely became stable dynamically (and it is still unstable thermodynamically) at 900, while it is becoming increasingly stable with temperature increasing to 1500 K. The effect can also be seen as a "precursor" showing that the temperature dependence of physical properties in bcc Ti might be anomalous. Let us demonstrate that this is indeed the case considering temperature dependence of elastic properties of β-Ti.

8
The calculated elastic constants and polycrystalline moduli are presented in Fig. 4. MTP predicts C ij s with negligible statistical error and in excellent agreement with results of AIMD simulations [35]. Note that because of the harder convergence criteria, the 220 static DFT calculations used for training MTP in terms of the computational time corresponds to about 1000 AIMD time steps. Thus, we estimate that our scheme allows for over two orders of magnitude speedup compared to calculating elastic constants from AIMD with the same accuracy and reliability of the calculated results. Analyzing the calculated elastic moduli, In experimentally synthesized GUM alloys the Elinvar behavior is explained through a strain glass state of the material with high defect concentration [31]. Random nanodomains of parent and martensite phases contribute opposite to elasticity, which results in that matreial's elastic moduli do not decrease with temperature. Interestingly, Elinvar effect in β-Ti predicted in our simulations is an intrinsic property of the material. Giving the fact that GUM alloys have complex chemical compositions the Elinvar property of elemental bcc-Ti observed in the wide temperature interval 1100-1700 K is unique.
In summary, we have shown that an efficient on-the-fly machine learning approach can be utilized to predict finite-temperature elastic behavior of materials with around two orders of