Moment tensor potential for static and dynamic investigations of screw dislocations in bcc Nb

A new machine-learning interatomic potential, specifically a moment tensor potential (MTP), is developed for the study of screw-dislocation properties in body-centered-cubic (bcc) Nb in the thermally- and stress-assisted temperature regime. Importantly, configurations with straight screw dislocations and with kink pairs are included in the training set. The resulting MTP reproduces with near density-functional theory (DFT) accuracy a broad range of physical properties of bcc Nb, in particular, the Peierls barrier and the compact screw-dislocation core structure. Moreover, it accurately reproduces the energy of the easy core and the twinning-anti-twinning asymmetry of the critical resolved shear stress (CRSS). Thereby, the developed MTP enables large-scale molecular dynamics simulations with near DFT accuracy of properties such as for example the Peierls stress, the critical waiting time for the onset of screw dislocation movement, atomic trajectories of screw dislocation migration, as well as the temperature dependence of the CRSS. A critical assessment of previous results obtained with classical embedded atom method potentials thus becomes possible.


Introduction
Plastic deformation in body-centered-cubic (bcc) metals is controlled mainly by the motion of screw dislocations.It is generally accepted that the screw dislocation motion in bcc metals is mediated by thermally-and stress-assisted kink-pair nucleation and growth [1].Detailed investigations of the screw dislocation mobility, the Gibbs energy and entropy of kinkpair formation were recently performed for the case of bcc Nb with standard and accelerated molecular dynamics (MD) utilizing an embedded-atom method (EAM) potential [2][3][4].Classical interatomic potentials, like EAM, are usually fitted directly to the bulk properties of the material and are computationally fast, but generally have limited capability of modelling defect properties because their analytical form is neither flexible nor easily extensible.For example for bcc Nb, the EAM potential of Farkas and Jones [5] predicts a too low Peierls barrier compared to previous ab initio calculations, and an enthalpy of kink pair formation which is an order of magnitude smaller than experimental values [2,3].
Ab initio methods, such as density-functional-theory (DFT), provide the necessary accuracy in the atomic interactions.However, dynamical simulations of screw dislocations require extended time scales and large supercells due to the long-range stress fields of the screw dislocations [1].Such simulations are still beyond the capability of current DFT methods, which are limited to a few hundred of atoms and several thousand of time steps (picoseconds) by their computational demands.Machine learning, on the other hand, has emerged as a promising new approach in material science for the construction of interatomic potentials, combining DFT accuracy with the computational speed comparable to that of classical empirical potentials.The progress in this fast-growing field has been reviewed in several recent articles [6][7][8][9][10].A key ingredient of all machine-learning interatomic potentials (MLIPs), regardless of their specific functional form, is the fitting (training) of the MLIP to a training set of energies, forces and/or stresses, generated using DFT or other quantum-mechanical (QM) methods [6][7][8].
Despite the great recent success, the development of general-purpose MLIPs remains a very challenging task [8].This is the reason why most MLIPs are developed as 'specific purpose' potentials and have poor transferability to atomic environments that are not included in the training set [9,10].This is also the case for bcc Nb [11][12][13][14][15].A spectral neighbour analysis potential (SNAP) [11] was developed for studying strengthening mechanisms in NbMoTaW multi-principal element alloys.A Gaussian approximation potential (GAP) for Nb was augmented with an analytical repulsive potential for radiation damage simulations involving highenergy collisions [12].A moment tensor potential (MTP) [13] was developed for the predication of the unstable stacking-fault energy in Nb-containing random alloys.Another MTP was developed in [14] for the description of hydrogen diffusivity in bcc metals.Yet other MTPs for bcc metals were developed recently to aid the prediction of accurate thermodynamic properties [15,16].The respective MTP for bcc Nb accurately predicts the vibrational part of the heat capacity and the bulk modulus at high temperatures up to the melting point [17], but the accuracy in reproducing the low temperature phonon dispersion is limited [15], which implies an inaccuracy in the elastic constants.
None of these MLIPs was specifically trained for the study of screw dislocations in bcc Nb and the corresponding training sets did not contain configurations along the Peierls barrier profile.The SNAP for example leads to a non-degenerate screw dislocation core [11].Apart from that, no data for the mobility of screw dislocations in bcc Nb using these MLIPs has been reported in the original papers.Therefore, in the present paper, we report on the development of a new machine-learning potential for bcc Nb, specifically trained for the study of screw dislocation mobility in the thermally-and stress-assisted temperature regime at low temperatures (⩽300 K).
Bcc Nb represents a very challenging test case for the construction of such an MLIP.Nb exhibits several elastic anomalies ( [2,18]): a low c 44 elastic constant, the lowest elastic anisotropy (Zener ratio), defined as A = 2c 44 /(c 11 − c 12 ), and the lowest G 111 shear modulus (G 111 = (c 11 − c 12 + c 44 )/3) among all bcc transition metals.Moreover, c 44 goes through a minimum at about 400-500 K, which is believed to be a Fermi surface effect [19].The phonon dispersion curves show some very unusual features, for example, the crossing of the longitudinal and the transverse branches near the H point [20,21], which are not observed for other bcc metals like W, Mo and Fe (similar features are also observed in Ta, but are not so pronounced as in Nb [22]).The interactions in bcc Nb are very long-ranged and the proper fitting of the phonon dispersion curves requires taking into account up to the 8th nearest-neighbour [21].Nb has also the lowest enthalpy of kink pair formation, compared to the other bcc transition metals [23].
To tackle these challenges inherent to bcc Nb, we use the MTP formalism developed by Shapeev [24].The fitting of the MTP for a single-component system (like bcc Nb) is fast and does not require much data, compared to other MLIPs, because the MTP utilizes a compact basis set with not too many parameters (compared to, e.g.GAP and artificial neural networks (ANN) [24]) and the evaluation cost is independent of the number of neighbours [9].In addition, MTP generally shows an excellent trade-off between execution speed and accuracy, compared to other MLIPs [25].

MTP
Machine-learning potentials have claimed their place as a more precise (though more costly) counterpart of their (semi-)empirical predecessors.Here, we rely on the MTPs as they are cost-effective compared to other MLIPs, require not so much data for training, compared to MLIPs like GAP and ANN, and since they were successfully applied to the study of metallic alloys [26][27][28][29].The theoretical description of the MTP formalism is given in [24].As common, the optimal MTP parameters are determined by minimizing the objective function, defined as the weighted sum of squared discrepancies between the DFT and MTP energies, forces and virial stresses for all atomic configurations in the training set.The weights of the energy, force and stress contributions to the objective function are set equal to 1.0, 0.01 Å 2 and 0.001 Å 6 , respectively.The optimal parameters are determined iteratively using the Broyden-Fletcher-Goldfarb-Shanno algorithm, starting from randomly initialized MTP parameters.The MLIP package [30] is used for the generation of the MTPs.Our results show that the relevant Nb properties can significantly vary among MTPs even if the potentials exhibit similar root-meansquare (RMS) energy errors on the training set.Potentials with smaller force errors, however, tend to provide better results.We have tested many different combinations of training configurations, cutoff radii, weighting schemes and MTPs of different levels, before choosing the final MTP of level 16 (see [24] for the level definition) with a cutoff radius of 5 Å (which is between the 3rd and the 4th nearest-neighbour in bcc Nb) and 125 fitted coefficients in total.

Training set
The assembly of an optimized training set is a key point for the development of an accurate and robust MLIP for a given application [6][7][8][9].MLIP fitting aims at reproducing the potential energy surface of the material, and not directly a given set of relevant physical properties.The relation between the RMS energy errors and the accuracy in predicting specific material properties is still an open question in the MLIP theory [9].In the case of MTPs additionally, it is not possible (like, e.g. for Gaussian processes [31]) to calculate explicitly the correlation between different configurations in the set.The question of a proper training set is thus always central to the construction of an MLIP, as it will fail to provide accurate predictions for underrepresented atomic environments.The so-called active learning approach can help to optimize the training set [26,32], however it was recently reported that active learning needs to be supplemented with human intuition to provide the best results [13,27].That is why in the present paper, based on previous experience with dislocation simulations [2][3][4], we use a large hand-made data set of reference structures judiciously optimized using physical insight and additionally by trial-and-error, targeted to describe well simultaneously the elastic, thermodynamic and defect properties of bcc Nb at low temperatures as well as the static screw-dislocation properties.Other MTPs for bcc Nb employed smaller training sets by using the active learning approach [13,14].A larger set takes more time to be created, but can be expected to lead to an MLIP with a broader applicability range.
The final composition of our training set, described below, is thus a result of a very thorough investigation.First of all, one cannot simply include realistic dislocation configurations, as they are not feasible for DFT.The task is to compose a proper training set, including relevant atomic environments, allowing for the MLIP to predict the correct forces in the dislocation region.Having a too simple training set with only few structure types (e.g.bcc, fcc, and a few deformed configurations) leads to a bad prediction of the dislocation properties, despite the fact that the training errors look very promising.On the other hand, taking too strongly diverse data leads to an MLIP with intolerably large training errors.In particular, inclusion of some properties that a priori may be considered useful (e.g.Bain deformation path at non-equilibrium volumes, tetragonal path plus random atomic displacements, ideal shear strength) leads to a significant increase of the RMS errors.Our final, optimized composition of the training set allows us to fit a relatively light-weight and fast, yet accurate MLIP (see sections 2.1 and 4), compared to other MLIPs.
Our training set consists of 18 subsets with in total 316 individual atomic configurations.The subsets are classified for clarity into 4 weakly-correlated groups (see table 1), which target the description of different groups of properties: the potential energy of several possible phases of Nb (bcc, fcc, body-centered tetragonal (bct), simple cubic (sc)), the elastic properties, the stacking-fault energy surfaces and defect properties (vacancies, surfaces, screw dislocations).For some configurations, we have not fitted the forces, because they were either zero by symmetry or could not be fitted reliably.The details of each group are described below.

Group A.
The training subsets in group A aim at an accurate reproduction of the general potential energy surface by including: the bcc phase at different volumes (training subset NA1; cf table 1), the higher-energy fcc phase (subset NA2), the bct and the sc structure as well as different distorted structures along the volume-preserving tetragonal (subset NA3), orthorhombic (subset NA4) and trigonal (subset NA5) deformation paths.These deformation modes are characterized by a single deformation parameter (denoted c/a ratio in the case of the tetragonal distortion and p for the orthorhombic as well as for the trigonal deformation paths).The paths connect the bcc with the fcc, sc and bct structures through a series of nonequilibrium configurations.A more detailed description of these deformation modes is given, e.g. in [33].The correct reproduction of the energy difference between the bcc and the fcc structures (∆E fcc−bcc ) is crucial, because recently it has been shown that the screw dislocation core structure is governed by ∆E fcc−bcc ([34] and references therein).The orthorhombic deformation path (subset NA4), passing through the bct structure, is particularly important for the proper description of the Peierls barrier, because the saddle-point configurations for shearing along {211}⟨111⟩ and {110}⟨111⟩ both have the same bct structure [18,35].The trigonal path (subset NA5) is also important for shearing simulations, because it represents a homogeneous deformation corresponding to the extension along the [111] direction, while keeping the atomic volume fixed.The configurations in training subset NA6 are generated by first running MD simulations under NPT conditions with the Nose-Hoover thermostat at 300 K using an initial low-level MTP and then performing DFT calculations (see section 2.3) at selected snapshots without structural relaxation.The selection of the MD snapshots is done with the active learning algorithm [26,32], by sampling up to 20 configurations at each step.

Group B.
The stress and dilatation fields of 1 2 ⟨111⟩ screw dislocations depend sensitively on the elastic constants [36].The nucleation energy of a kink-pair is proportional to G 111 b 2 /2π [2,37], where b is the length of the Burgers vector.It is generally accepted that the properties of screw dislocations in bcc metals are closely related to and depend on the characteristics of the stacking-fault energy (SFE) surfaces [39,40].Correspondingly, group C of the training set targets the accurate reproduction of the (110) and the (211) stacking-fault energy curves (training subsets NC1 and NC2, respectively) along the [111] direction.Such configurations were also used previously for fitting the GAP [12] and the MTP [13] for bcc Nb.

Group D.
Group D targets the accurate description of various defect structures in bcc Nb that are of key relevance to the simulation of screw dislocations.It is generally accepted that the climb of screw dislocations is mediated by the presence of vacancies [1].Training subset ND1 contains, correspondingly, several configurations with an unrelaxed single vacancy.
Single-dislocation models with free boundary conditions, perpendicular to the dislocation line, are often used in MD simulations of screw dislocation glide.In order to account properly for the presence of free surfaces, group D contains atomic configurations (training subset ND2) with (100), (110), (111) or (211) free surfaces perpendicular to the Z direction and small random displacements of the atoms.For the DFT calculations using the ND2 configurations, a vacuum layer of 20 Å is added along the Z direction, perpendicular to the slabs.
Most importantly, we include in group D diverse models of different sizes containing screw dislocations, because inclusion in the training set of only the SFE curves is not sufficient to capture the strongly-anisotropic dislocation core properties.For the dipoles in subset ND3 we use an orthogonal simulation box with 4b length of the screw dislocation line.ND4 represents sheared atomic configurations containing straight dislocations and dislocations with single or multiple kink pairs (see figure 1).The length of the dislocation line is 3b.Vacuum layers of 20 Å are added along the X and Z directions for the DFT calculations.The atomic configurations in ND5 represent dipoles with 1b length of the dislocation line in a special quadrupole arrangement (see section 2.4 below) for which the elastic interactions between the two dislocations cancel out [41].Such a geometry has been commonly used in previous DFT calculations of the static screw dislocation properties of Nb and other bcc metals [42][43][44].The advantage of such configurations for DFT calculations is the small size and the tri-periodic boundary conditions.The ND5 configurations are generated by the nudged elastic band (NEB) method [45] along the minimum energy path between two easy-core dislocation positions using an initial set of replicas, constructed by linear interpolation.The ND5 training subset is specially included, because many empirical potentials for bcc metals often predict a 'double-hump' Peierls barrier, corresponding to a metastable core structure (see [42] and references therein).
Several MLIPs for other bcc metals (e.g.Fe, W and Ta) were fitted using training sets containing configurations with straight screw dislocations [47][48][49].Screw dislocations with kink pairs were not included into the corresponding training sets.To the best of our knowledge, this is the first machine-learning potential for bcc Nb, the training set of which contains screwdislocation configurations with and without kink pairs.

DFT calculations
The DFT calculations were performed with the Vienna Ab initio Simulation Package (VASP [50]) using the Nb_pv projector-augmented wave (PAW) potential for Nb with the valence configuration of 4p 6 4d 4 5 s 1 (11 valence electrons), distributed with VASP.Except for the energy of the isolated Nb atom, we used non-spin-polarized calculations.We employed the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) exchangecorrelation functional.On the basis of convergence tests for the cohesive energy, the planewave expansion cutoff energy was set to 520 eV for all DFT calculations.Such a high cutoff energy ensures a high accuracy for the components of the stress tensor as well.To sample the Brillouin zone, we used the Monkhorst-Pack method [51].For groups A-C as well as the ND1 training subset, we used a k-point grid spacing of 0.075 Å −1 .For the larger configurations in training subsets ND2-ND5 we used smaller k-points grids due to computational costs.The first-order Methfessel-Paxton method [52] with a smearing width of 0.2 eV was used for integration in the Brillouin zone.During electronic relaxation, convergence was assumed when the energy variation between two subsequent electronic self-consistent steps was below 10 −4 eV.

Static and MD simulations
The MTP-based static and MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [53], compiled together with the MLIP module.Most of the screw-dislocation models were generated using the program ATOMSK [46].An in-house program was used for the generation of a periodic array of dislocations (PADs).Atomic structure optimization was performed with the conjugate gradient method.Convergence in LAMMPS was assumed when the energy differences dropped below 10 −8 eV and the force differences below 10 −8 eV Å −1 .During the MD simulations at finite temperatures, the system was initially equillibrated at the desired temperature using the Nose-Hoover thermostat, starting from an ensemble of velocities with a uniform distribution under a constant number, volume and temperature integration (NVT conditions).The time step for the MD simulations was 1 fs.
For the dislocation models with (110) maximum resolved shear stress (MRSS) plane, the orientation of the supercell was X || For the static calculations of the Peierls stress using large supercells, energy and force minimization was first performed in the unstressed state, using the conjugate gradient algorithm at 0 K.The shear strain along the [111] direction was then increased incrementally and an energy/force minimization was carried out at every step with LAMMPS [53].For every resulting configuration, the position of the screw dislocation on the (111) plane was determined using the dislocation analysis, as implemented in the OVITO program [93].The shear strain and the corresponding shear stress, at which the screw dislocation started to move, were taken as the critical strain and the critical stress (Peierls stress).Large conventional supercells, containing from 686 to 85 750 atoms, were used to check the convergence and the effect of vacancy-vacancy image interactions on both the unrelaxed and relaxed single-vacancy formation energy (E vac ).For the calculation of the phonon dispersion curves with the MTP, we used MD simulations and a 12 × 12 × 12 primitive supercell, containing 1728 atoms; 1 × 10 6 time steps for equilibration and 3 × 10 6 times steps for the determination of the phonon frequencies.The dynamical matrix was obtained from the thermal Green's function according to the fluctuation-dissipation theorem as implemented in LAMMPS [54,55].The phonon dispersion curves were calculated from the dynamical matrix using the program PHANA [56].The mean square displacement ⟨u 2 ⟩ at different temperatures was calculated using MD simulations (1 × 10 4 time steps for equilibration and 5 × 10 4 times steps for the averaging) in a 9 × 9 × 9 conventional supercell.We used the following thermo-elastic expression for the calculation of the thermal expansion coefficient α at constant volume: where B is the bulk modulus and p is the pressure.MD simulations using a 9 × 9 × 9 conventional supercell were performed at several temperatures below and above 300 K. We found a linear temperature dependence of the pressure, from which α was determined using equation (1).

Validation of the MTP
An accurate and unbiased agreement between the DFT and MTP energies and forces is demonstrated in figure 2. The energy RMS error for the selected MTP is 4.38 meV (figure 2(a)) and the force-component RMS error is 59 meV Å −1 (figure 2(b)), the RMS error for the stresses is 8.7 meV Å −3 .These RMS errors are of the order of most MLIPs [8].Taking into account the diversity of the training set, this is a decent result for an MLIP with only 125 coefficients and 90 µs/atom inference time on a single core (see section 3.3).A broad range of lattice and thermodynamic properties was calculated for further validation of the developed MTP.Table 2 summarizes the basic bulk properties of bcc Nb.The present MTP results are compared with our DFT calculations, included in the training set (groups A  [13]. c [11].d [12].e [57].f [58].g [59].h [60].i [61], at 4.2 K. j [62].k Calculated in the present study using the MTP of [13].l Calculated as: and B, see table 1), other MLIPs, previous electronic-structure calculations and experimental data.Our MTP reproduces well the DFT equilibrium lattice constant, the cohesive energies E coh of the bcc and fcc phases (see figure 3) as well as the energy difference between the fcc and the bcc structures ∆E fcc−bcc .The elastic constants are also in good agreement with our and previous [58] DFT calculations.Similar results for c 11 and c 12 are obtained using GAP [12], but GAP [12] significantly underestimates c 44 .Both the MTP, developed in [14], and the SNAP from [11] overestimate c 11 , while the MTP, from [13], underestimates c 11 and overestimates c 12 as well as the Zener ratio.All MLIPs for bcc Nb underestimate the c 44 constant, which  reflects the inability of the underlying DFT methods to correctly predict c 44 , regardless of the exchange-correlation functional and the approach to solve the Kohn-Sham equations.
The energies along the tetragonal, orthorhombic and trigonal deformation paths (included in group A, see table 1), calculated with the MTP, are in excellent agreement with our DFT calculations (figure 4) even at large deformations.Figure 4 indicates that the MTP accurately samples the distorted atomic configurations between the simple bcc, fcc, bct and sc structures.
The SFE curves, calculated with our MTP, are also in good agreement with the DFT calculations (see figure 5).The (110) MTP curve is symmetric, while the (211) MTP curve is asymmetric with respect to the middle point (displacement 0.5b), as required by symmetry.These calculations are based on 8 atomic layers, perpendicular to the slip plane.Extensive additional calculations using larger stacks of up to 24 atomic layers indicate that the SFE curves remain practically the same.No local minima, characteristic for metastable stacking faults, are observed.The energies of the unstable stacking fault γ us , defined as the SFE maximum at a displacement 0.5b, are compared in table 3 with previous results.The unrelaxed   [11]. 3 [13]. 4 [64]. 5 [58].
SFEs reported in [63] are slightly higher than our DFT results, probably because these authors used the LDA approximation.Although the relaxed SFEs were not included in the training set, our MTP results are in relatively good agreement with previous relaxed DFT calculations [58], demonstrating that the training set contains relevant information (environments) for the MTP to carry out the desired simulations without (or just with little) extrapolation.On the other hand, the MTP of [13] underestimates, while the SNAP [11], used in [64], significantly overestimates the relaxed stacking fault energies.
In order to examine our MTP at finite temperatures we have performed extensive MD simulations at 300 K.The MTP phonon dispersion curves at 300 K along the main symmetry directions are compared to experiments [20] in figure 6.The MTP reproduces correctly the shape and all anomalies in the dispersion curves of bcc Nb, although the frequencies of some zone-boundary phonons deviate slightly from the experimental results.Several other vibrational properties of bcc Nb, not represented in the training set, are compared with experiments in table 4. The Grüneisen parameter is calculated as αV mol B/C p , where V mol is the molar volume and C p is the heat capacity at constant pressure [65].In table 4 our MTP data are closer to experiment than the GAP and tight binding counterparts.The bcc structure of Nb remains stable using the developed MTP up to 2000 K at ambient pressure and up to 50 GPa at   [65], at 1500 K. c [68].d [69].300 K, although our MTP was not specifically trained at such high temperatures and pressures.This demonstrates further the robustness of the MTP.
The relative errors between the MTP and DFT are summarised in figure 7. The developed MTP achieves DFT accuracy of about ±10% for properties included in the training set (figure 7(a)), and errors of ±15% for properties not directly included in the training set (figure 7(b)).Despite that the trend (larger error for the out-of-training properties) is natural, the individual error distribution for the different properties is hard to predict; i.e. some of the properties on which we trained are reproduced worse than some of the properties not included in the training set (e.g.c 44 versus c 12 ).

Results and discussion
In this section, we provide a comprehensive study of the screw-dislocation properties in bcc Nb using large supercells, taking advantage of the computational efficiency of the generated and thoroughly validated MTP (section 2.5).Sections 3.1 and 3.2 contain static and dynamic screwdislocation properties, respectively, while section 3.3 presents CPU timings for the large-scale dislocation simulations with the MTP.

Static screw-dislocation properties
The core structure of the 1  2 ⟨111⟩ screw dislocation in bcc Nb, relaxed with our MTP at 0 K and zero applied shear stress, is shown in figure 8 as a differential displacement (DD) map [70] obtained with the program DDPLOT [71].The DD plot shows a non-planar core structure without any indication for splitting along any of the three [211] directions and thus corresponds to the so-called non-degenerate (compact) type of easy core.This is in qualitative agreement with previous calculations using DFT [42][43][44]72] and empirical potentials [2,40,72,73] for bcc metals.
The energy of the easy core E core was determined using (110) dipole models of increasing size along X and the equation [72]: where ∆E is the energy difference (per dislocation and per unit dislocation length) between the total energies of the dipole model and the corresponding perfect lattice, K s is a shear modulus, depending on the modified elastic compliances [2], d is the distance between the two dislocations along the X axis, r core is the inner radius of the core.Due to the long-range nature of the interactions in bcc Nb, a core radius equal to the 3rd nearest-neighbour distance (r core ∼ 1.633b) was used in equation (2).A linear dependence of ∆E versus ln(d/r core ) was observed for several dipole models with d ranging from 39 to 161 Å.The easy-core energy, determined from the intercept of a linear fit using equation ( 2), is equal to 0.22 (12) eV/b and the shear modulus, determined from the slope, is equal to 30.5(1.0)GPa.The MTP easy-core energy is in very good agreement with previous DFT calculations (see table 5).The value K s = 30.5(1.0)GPa is in good agreement with the G ′ 4 shear modulus, calculated directly from the elastic constants (see table 2).
The ability of the screw dislocations to glide on the (110) plane is governed by the energy barrier between two adjacent easy-core positions, known as Peierls barrier (Peierls potential).The Peierls barrier was computed with the MTP using a skewed supercell containing 231 Figure 8. Differential displacement map of the 1  2 ⟨111⟩ screw dislocation obtained with the present MTP, projected on the (111) plane using the DDPLOT program [71].The length of the arrows in the plot is proportional to the atomic displacement of the atoms in the core, relative to their positions in the ideal crystal.The atoms belonging to different layers in the ideal crystal, perpendicular to the [111] zone axis, are shown in different grey colours.atoms (see section 2.4) and the nudged elastic band (NEB) method [45], as implemented in LAMMPS [53].The calculations were performed with the FIRE minimization algorithm and a parallel NEB spring constant of 5.0 eV Å −2 .Additional tests with different parallel spring constants from 1.0 to 10.0 eV Å −2 showed that the choice of the spring constant does not affect the NEB results.The MTP and our DFT calculations give similar results for the energy along the transition path, showing a single maximum (figure 9).The height of the barrier is in close agreement with previous DFT calculations (see table 5).The EAM potential [5] also leads to a 'single-hump' profile, but the Peierls barrier is lower (∼19 meV/b) [2].On the other hand, NEB calculations, performed in the present study, using the MTP from [13] and the SNAP of [11] both lead to a 'double-hump' Peierls profile.The Peierls stress was calculated with the MTP for all three slip planes, (110), ( 211) and (123), using a range of supercells of different sizes, shapes and boundary conditions.Previous DFT calculations for bcc Nb were done only for the (110) slip plane and very small skewed supercells (135 atoms in [43] and 231 atoms in [42]).In order to compare with the results published in [42], we used the same skewed supercell with quadrupole setting [41] (see section 2.5).Stress was applied by shearing the supercell, i.e. by increasing the a 2 basis vector, a [41,74].For every x value, the atoms in the supercell were relaxed.A linear increase of the shear stress was observed with increasing x up to a critical value x crit (similar to bcc Fe [74]), at which the two dislocations simultaneously jumped to the next Peierls valley.The Peierls stress was taken equal to the shear stress at x crit .The (110) Peierls stress, determined with the MTP in this way, is higher than our DFT value as well as the Peierls stress reported in [42] (see table 6).
Wang et al [64] reported Peierls stresses for the (110), ( 211) and (123) slip planes using the SNAP, developed in [11], and PADs containing about 50 000 atoms.We determined the Figure 10.Dependence of the (211) Peierls stresses for large atomic configurations, containing a single dislocation (SD), on the length of the supercell along the X axis (Lx): full circles-twining orientation with positive sense of shearing (T+), full squaresanti-twining orientation with negative sense of shearing (AT−), empty circles-twining orientation with negative sense of shearing (T−), empty triangles-anti-twinning orientation with positive sense of shearing (AT+).The lines are a guide for the eye.
Peierls stresses using PADs of the same size, as in [64], for comparison.The SNAP Peierls stresses for the (211) plane are similar to our MTP results, but the SNAP Peierls stresses for the (110) and the (123) slip planes deviate significantly from our MTP values (table 6).We tentatively relate these differences to the fact that the SNAP [11], used in [64], overestimates the c 11 and c 12 elastic constants (see table 2) and overestimates the unstable SFEs for all three slip planes (see table 3).
Three types of boundary conditions are commonly used in the simulations of screw dislocation mobility: (a) dipoles with periodic boundary conditions along all three directions; (b) PADs with periodic boundary conditions along the glide direction and the dislocation line as well as free boundary conditions perpendicular to the glide plane; (c) single dislocation (SD) models with periodic boundary conditions only along the dislocation line.
We have performed extensive simulations using large supercells (compared to previous DFT calculations [42][43][44]) to investigate the effect of supercell size along the glide direction (L x ) on the Peierls stress for different boundary conditions.The Peierls stresses for the (110) slip plane with positive and negative sense of shearing, obtained with our MTP, are the same, as required by symmetry for all supercell sizes (see table 6) and different boundary conditions.The MTP also correctly reproduces the intrinsic twinning-anti-twinning asymmetry of the Peierls stresses for the (211) slip plane and different senses of shearing [40].
The Peierls stresses for the (110) slip plane are well converged already for L x ⩾ 100 Å, regardless of the boundary conditions used.For the (211) slip plane and the SD models, the Peierls stresses are well converged for L x ⩾ 280 Å (see figure 10). Figure 10 also demonstrates that the Peierls stresses, calculated with our MTP, comply very well with the important symmetry requirement [40] that the stresses for the twinning configuration and a positive sense of shearing should be equal, within numerical error limits, to the stresses for the anti-twinning configuration and a negative sense of shearing.For the (211) dipoles and PADs, convergence is also observed for L x ⩾ 280 Å.We have calculated the (110) Peierls stress using PADs with different dislocation-line lengths from 1b to 80b.No dependence of the (110) Peierls stress on the dislocation line length is observed with the MTP, contrary to the results obtained with SNAP [64].
The average Peierls stress for the (110) slip plane is practically independent of the boundary conditions within statistical limits, while a small increase is observed for the (211) Peierls stress when going from dipoles, to PADs and SDs.Although the unstable stacking-fault energy for the (211) plane is larger than that for the (110) plane (see table 3), there is a correlation of the SFE with the Peierls stresses only for the SD configurations.A similar observation was reported for other bcc metals as well [75].
The twinning-anti-twining asymmetry of the {211} Peierls stresses in bcc metals is also manifested in the fact that the Peierls stresses do not obey the Schmid law, according to which the CRSS should be inversely proportional to cos(χ) (dashed line in figure 11), where χ is the angle between a given slip plane and the (10 1) slip plane.The value χ = −30 • corresponds to the ( 211) twinning slip plane and χ = 30 • to the anti-twinning ( 11 2) slip plane.Using PAD models with about 50 000-54 000 atoms and special crystallographic orientations along the normal to the glide plane and along the glide direction, we have determined the Peierls stress as a function of the χ angle (see figure 11).The orientations of the corresponding supercells are given in the appendix.The atomistic simulations, visualized with OVITO [93], show that for all values χ < 26 • , the screw dislocation starts to move parallel to the ( 110) plane.Correspondingly, the CRSS values can be well fitted (full line in figure 11) with the phenomenological equation, proposed by Gröger et al [76] where τ * is the non-glide stress acting perpendicular to the screw dislocation, τ c is the CRSS at χ = 0 • and a 1 , a 2 and a 3 are fit parameters.The predicted CRSS-χ dependence clearly deviates from the Schimd law.Table 5 and figure 11 demonstrate that most of the static screw-dislocation properties, calculated with our MTP, comply very well with the symmetry of the bcc lattice and are in good agreement with our and previous DFT calculations.In particular, the MTP predicts simultaneously a non-degenerate dislocation core and an accurate Peierls potential, which is not the case for many other interatomic potentials including MLIPs for other bcc metals [34].However, the Peierls stresses for both the (110) and the (211) slip planes are about 3-4 times larger than the experimentally reported Peierls stress of about 0.4 GPa [77].A significant overestimation of the Peierls stresses, compared to experimental values, is also reported for various other bcc metals using SNAP [11,64], for bcc V using both GAP [12] and a deep MLIP framework [34] as well as for bcc Nb using an angular dependent potential (ADP) [78].These observations show that DFT accuracy (±10%-15%), achieved for most of the properties of bcc metals using a carefully-trained MLIP, does not lead automatically to Peierls stresses comparable to experiments.Further studies are required to clarify the discrepancy.

Dynamic screw-dislocation properties
Very few data exist for the glide motion of screw dislocations in pure bcc metals using large supercells and MLIPs [31,34,79].Both Wang et al [34], using an extended modified EAM potential for bcc V at 77 K, and Maresca et al [79], using a GAP potential for bcc Fe at 200 K, reported that the screw dislocation starts to glide on the (101) plane via double-kink nucleation and growth.However, no MLIP data are available for bcc Nb regarding, in particular, screw dislocation trajectories over extended time scales and the dependence of the CRSS on strain rate and temperature.To simulate screw dislocation glide we used PADs, because this setup allows to analyze the dislocation trajectory over extended MD time scales, keeping the size of the supercell along the glide direction relatively small.A similar approach was applied to bcc Fe using a GAP [79].Different computational schemes (constant-stress or constant-strain loading) can be used for shearing the screw dislocation models.In experiments, tensile or shear deformation is performed at constant strain rate [80].Moreover, our own MD simulations for the screw dislocation in bcc Nb and published work for an edge dislocation in fcc Cu [81] have shown that the constant-stress loading mode leads initially to shear-stress oscillations using PADs.For these reasons, we used PADs (containing 96 000 atoms) and the constantstrain loading mode under NVE conditions.In these simulations we used a positive sense of shearing.Multiple runs with different random numbers during equilibration were performed at selected temperatures for better statistics.The temperature of the mobile atoms remained constant within ±1 K under these conditions for all strain rates, demonstrating further the robustness of the developed MTP.

Strain rate effects.
The extensive literature on the mechanical properties of bcc Nb, studied experimentally at different strain rates up to 1.1 × 10 3 s −1 , was recently reviewed in [82].In particular, Seeger and Holzwart [80] reported data for the strain-rate dependence of the CRSS down to 120 K. Their results show that in the low strain-rate regime (rates between 6.5 × 10 −5 and 3.5 × 10 −3 s −1 ) the CRSS (flow stress) of bcc Nb increases with increasing strain rate.The strain rates in MD simulations have to be much higher in order to achieve meaningful shear strains in reasonable computational time.The developed MTP was tested and shows good robustness for strain rates up to 5 × 10 9 s −1 .In this high strain-rate regime, our MD simulations at 5 K show that the CRSS (τ c ) is different for the two slip planes, 1.38 ± 0.1 GPa for the (211) PAD models and 1.55 ± 0.05 GPa for the (110) PAD models, but practically constant with increasing strain rate suggesting that different mechanisms of screw dislocation mobility might be operational at low and high strain rates.The CRSS τ c is commonly expressed as where γ c is the critical strain at which the screw dislocation starts to move and G is the shear modulus.From equation (4), the critical waiting time t c for the onset of screw dislocation movement should decrease inversely with increasing strain rate at constant CRSS (taking into account that γ c = γt c ) where t 0 is a fit parameter and γ is the strain rate.Indeed, plotting t c versus the strain rate on a log-log scale, a linear increase of t c with decreasing strain rate is observed for both the (211) and the (110) models (see figure 12).The waiting times for the (110) model are throughout slightly larger than for the (211) model (at the same strain rate), but the strain-rate dependence of the critical waiting time (slopes in figure 12) is practically independent of the slip plane.

Trajectories and slip planes.
Identification of the most probable slip plane in bcc metals is still a matter of debate (see for example [75,80,83]).As for bcc Nb, contradictory data have been published regarding the favorable slip plane.Plastic deformation was reported by many authors to occur at low temperatures predominantly on the (110) slip plane [83].But several authors [84,85] reported slip in bcc Nb on either the (110) or the (211) slip plane, depending on temperature and/or loading direction.Seeger and Holzwart [80] performed a detailed analysis of the temperature dependence of the CRSS (in the range 125-350 K) using different models of kink-pair nucleation, and concluded that their data is in quantitative agreement only with (211) as the elementary slip plane.MD simulations for bcc Nb [2], using the EAM potential of Farkas and Jones [5] and SD as well as PAD models, showed that the screw dislocation trajectories depend on the MRSS plane and that the screw dislocation moves by alternating slip on different (110) slip planes, leading to an effective (211) glide plane, independent of temperature.
The trajectories of the screw dislocation, obtained from the present MTP MD simulations, show already at 5 K a complex behaviour.The screw dislocation trajectories, when resolved in atomistic detail, are not parallel to the corresponding MRSS plane (see table 7 and figure 13).
In the case of the (110) MRSS plane, the effective glide direction is initially parallel to the (0 11) slip plane at about −60 • to the MRSS plane (see figure 13(a)) for all strain rates.After multiple cross-slip events, the screw dislocation starts to move effectively parallel to the (1 10) slip plane in the subsequent traverses of the supercell, but the trajectory is rather wavy (figure 13(b)).With decreasing strain rate, the cross-over from an effective glide along the (0 11) slip plane to an effective glide along the (1 10) slip plane takes place later and the waviness of the trajectory increases.On the contrary, for the Farkas and Jones EAM potential, the initial effective trajectory is parallel to the (1 21) slip plane at about −30 • to the MRSS plane [2].
In the case of the (211) MRSS plane, the effective glide direction is initially parallel to the ( 110) slip plane, at about −30 • to the MRSS plane (figure 13(c)) for all the studied strain rates.After multiple cross-slip events, the screw dislocation starts to move effectively parallel to the (2 11 ) slip plane in the subsequent traverses of the supercell (figure 13(d)).With decreasing strain rate, the cross-over from an effective glide along the ( 110) slip plane to an effective glide along the (2 11 ) slip plane takes place later.On the contrary, for the Farkas and Jones EAM potential, the screw dislocation starts to move immediately parallel to the (211) MRSS plane [2].
The screw dislocation trajectories also change with temperature.Above 5 K, the initial effective glide parallel to the ( 110) plane for the (211) PADs persists only for a short time due to the enhanced thermal fluctuations.After that the screw dislocation starts to glide parallel to the (2 11 ) plane, but the elementary jumps are not well defined and the trajectories are rather wavy (see figure 14).A similar trend with increasing temperature is also observed in the case of the (110) PADs.
The effective trajectories in the MD simulations consist of elementary jumps, parallel to specific planes (summarized in table 7).Initially, the MTP elementary jumps are predominantly parallel to alternating {211} symmetry-equivalent planes: (1 21)+( 112) in the case of the (1 10) models (figure 13  Using the present MD simulations with a DFT-accurate MLIP, it seems thus possible to reconcile the different observations of slip in bcc Nb at low temperatures.The effective glide (effective trajectory) is either on the (110) or (211) slip plane, depending on the loading direction, in agreement with the macroscopic slip-trace observations performed in the 1960s, but the elementary jumps are initially along the {211} slip planes, in agreement with the detailed analysis of Seeger and Holzwart [80].

3.2.3.
Temperature dependence of the CRSS.Already at 5 K the onset of the screw dislocation movement takes place by the nucleation and growth of a single kink pair (see figure 15), in agreement with previous theories [23] and similar to the results for bcc Nb [2] and other bcc metals (e.g.bcc Fe [86][87][88]) obtained with EAM potentials.To assess the thermally-activated regime of kink-pair nucleation and growth with the developed MTP over a broader temperature range, the CRSS was determined as a function of temperature, utilizing a strain rate of 6.18 × 10 8 s −1 and the (211) PAD model.The nucleation and growth of a critical (stable) kink pair at the onset of the screw dislocation glide is clearly observed up to 250-300 K with an average length of the critical kink pair equal to 5 ± 3 Å.At still higher temperatures the screw dislocation line becomes rather wavy and nucleation of multiple kink pairs takes place.The CRSS as a function of temperature is compared in figure 16 to literature data [2] obtained with the EAM potential of Farkas and Jones [5].The CRSS decreases initially with increasing temperature, forms a small 'hump' in the temperature range 250-350 K and reaches a constant athermal stress level τ a of about 0.55 GPa for T ⩾ 400 K.The presence of a 'hump' in  the CRSS-temperature relationship is predicted by the theory of Seeger [23].In the thermallyactivated regime, the temperature dependence of the CRSS can be phenomenologically written in the form [2,89] where ∆H kp (0) is the enthalpy of kink-pair formation at zero shear stress, τ 0 , p and q are phenomenological model parameters.The constant strain-rate factor γ0 is proportional to the dislocation density, the glide area per activation event and the attempt frequency [89], and it was taken equal to γ0 = 1.23 × 10 9 s −1 in the present case.CRSS data in figure 16, obtained with the present MTP, can be approximately modelled with equation ( 6) with parameters: τ a = 0.55 ± 0.02 GPa,τ 0 = 0.90 ± 0.02 GPa, ∆H kp (0) = 0.027 ± 0.003 eV, p = 0.5 and q = 1.25.This demonstrates further that the MTP leads to thermally-activated kink-pair nucleation and growth.
The CRSS curve, determined with the EAM potential, has overall a similar temperature dependence, but the 'hump' is not so well expressed and all CRSS values are about 50% lower.The corresponding model parameters are: τ a = 0.22 0.05 GPa, τ 0 = 0.75 GPa, ∆H kp (0) = 0.03 ± 0.01 eV, p = 0.5 and q = 1.45 [2].differences the CRSS values, obtained with the MTP and EAM potentials, can be elucidated taking into account the differences in the elementary jumps.Initially, the MTP elementary jumps are predominantly along the (211) slip planes (see above), while the EAM elementary jumps are predominantly along the (110) slip planes [2,3].The properties of the screw dislocations are closely related to the SFE curves [39,40], which characterize the resistance of a crystal to shear along a given crystallographic (slip) plane in the [111] direction.The MTP (211) SFE energy barrier is higher than the (211) SFE energy barrier obtained with the EAM potential [2].This consideration qualitatively explains why the CRSS values, obtained with the MTP, are higher than the CRSS values, obtained with the EAM potential (figure 16).
Most notably, however, the enthalpies of kink-pair nucleation, obtained with the MTP and EAM potentials, are very similar (∼0.03 eV) and much lower than the experimental value of 0.68 eV, reported by Seeger and Holzwart [80].This comparison indicates first of all that the low enthalpy of kink-pair nucleation, obtained earlier with the EAM potential, is not necessarily due to limitations of the EAM, as previously suggested [2].During the nucleation of the very first kink pair at 5 K with the MTP, the two individual kinks are always along {110} symmetry-equivalent slip planes, independent of the strain rate: along (1 10) + (0 11) slip planes in the case of (211) PAD models and along (0 11) + ( 101) slip planes in the case of (110) PAD models.The individual kinks during the nucleation of the very first kink pair are also along {110} slip planes when using the EAM potential [3].The same slip plane of the initial kinks qualitatively explains why the enthalpies of kink-pair formation, ∆H kp , are similar for the two potentials, despite their different accuracy.
As indicated in figure 15, the formation and evolution of the double kink is a complex cooperative event involving many atoms in the dislocation core.A direct computation of realistic kink-pair models is still beyond the capabilities of ab initio simulations.To circumvent this difficulty, Dezerald et al [44] employed small supercells with only 1 to 2 b length of the screw dislocation in order to derive DFT values for the line tension and the Peierls barrier, using NEB calculations at 0 K.These parameters were then used to calculate numerically the enthalpy of kink pair formation using a linearized version of the line-tension functional.Specifically, for bcc Nb they reported a value (1.28 eV), which is almost twice larger than the experimental value (0.64 eV [80]).We expect the here presented dynamic calculations at the relevant temperatures with a verified MLIP to better reflect experimental conditions.
The remaining discrepancy between experiment and simulations could originate in the different probed strain-rate regimes, manifested in different logarithmic strain-rate ratios − ln( γ γ0 ) (see equation ( 6)) in simulations (0.69 in the present case) and in experiments (30 for bcc Nb [90] and 31-32 for bcc Mo and W [76]).The recently proposed strain-reduction method combined with hyperdynamics [4] could provide a way to mitigate the strain-rate impact.The strongly differing dislocation densities between experiment and simulations could likewise contribute to the discrepancy.In any case, further investigations are required to elucidate the discrepancy.

Computational speed
We compare here the computational performance of the developed MTP with classical potentials and other MLIPs.Specifically, we compare with the EAM potential used in our previous MD simulations [2][3][4], the modified EAM potential developed by Lee and Baskes [91], the SNAP [11] and the GAP [12].In order to estimate the computational speed, we measured the time elapsed for running 10 000 MD steps on a 3 GHz Intel i5 processor with 6 cores under NPT conditions at 300 K using a 9 × 9 × 9 conventional supercell with 1458 atoms.The data in figure 17 shows that the speed of the MTP (in ns/day, as reported by LAMMPS) is 44 times larger than GAP, comparable to that of the SNAP, ∼12 times smaller than that of the MEAM potential and almost ∼37 times smaller than that of the EAM potential.We note that the current implementation of the MTP in the MLIP package is available on CPUs only.An implementation on graphical processing units (GPUs) could increase the performance of the MTP.

Summary and conclusions
A MTP has been developed for the accurate description of the complex behavior of bcc Nb in the low-temperature regime.A comprehensive validation of the MTP has been performed, showing that it reaches DFT accuracy to within ±10%-15% for a broad range of mechanical, thermodynamic and defect properties of bcc Nb.The MTP is especially suited for the modelling of screw dislocations, a property that has been achieved by including into the training set a diverse collection of screw-dislocation configurations, some of which contain kink pairs.The MTP can be used in the framework of the LAMMPS [53] code in a straightforward manner.
The MTP predicts correctly the non-degenerate screw dislocation core structure and reproduces very well the DFT energies of the easy-core configuration and of the Peierls barrier.The MTP also reproduces correctly the twinning-anti-twining asymmetry of the Peierls stress for the (211) slip plane as well as the orientation dependence of the CRSS at different values of the orientation angle χ between the ( 101) slip plane and the maximum resolved shear stress plane.The MTP thus enables large-scale atomistic simulations of screw dislocations at finite temperatures which are not feasible with first-principle methods.
The MTP demonstrates good robustness and extrapolation ability, preserving the stability of the bcc structure of Nb up to 2000 K at ambient pressure and up to 50 GPa pressure at 300 K, although such configurations were not included in the training set.We attribute this MTP robustness to the following two facts.Firstly, the functional form of our MTP (level 16) is relatively simple compared to other ML models, facilitating successful extrapolation to a good extent.Secondly, our hand-made training set contains enough diverse atomic structures (and, hence, atomic environments) of different origin, which provides a comprehensive sampling of the relevant part of the phase space.In combination, these two factors provide the MTP with the ability to model configurations which are not explicitly included in the training set.
In this respect, it is important to stress that we have deliberately avoided the active learning approach in the present research, instead relying on a careful by-hand design of the training set.This has given us full control over the relevant properties.In contrast, with active learning it is usually not clear what exactly is in the training set, and how it affects the relevant properties.Moreover, for better results active learning should be anyway reinforced with physical knowledge [13,27].
Compared to the other three previously published Nb MTPs, we used a more dedicated training set which results in a better reproduction of most properties of bcc Nb at low temperatures.The GAP [12] and its faster counterpart SNAP [11] also demonstrate good correlation with some of the DFT and experimental data for bcc Nb.Their less accurate prediction of other properties (especially the higher stacking fault energies and the 'double-hump' Peierls barrier), compared to the present MTP, can be attributed to our more precisely tailored training set.The slower computational speed of GAP [25] and the speed drop of SNAP with increasing number of components [92] (unlike MTP) render the MTP more appealing for simulations with several components and large supercells ( ≳100 000 atoms).The application of the developed MTP to large supercells has been already demonstrated in the present study.
The Peierls stresses for all three slip systems, (110), ( 211) and ( 123), obtained with the MTP, are about 3 times larger than the experimentally reported value for bcc Nb [77,80].This emphasises that the DFT accuracy of an MLIP does not necessarily ensure good agreement with experiment for all screw dislocation properties.
The CRSS remains practically constant for strain rates from 6.2 × 10 7 to 2.5 × 10 9 s −1 , leading to an increase of the critical waiting time for the onset of the screw dislocation movement with decreasing strain rate.The initial effective screw-dislocation trajectories for both the (110) and (211) MRSS planes are parallel to one of the (110) symmetry-equivalent slip planes, but the elementary slip jumps are predominantly along the (211) slip planes.This finding is in contrast to the results obtained earlier with a classical EAM potential [2], which had not been specially trained to reproduce accurately the SFE curves and the Peierls barrier.
The MTP allowed us to reveal that the time of the cross-over from the initial effective (110) MRSS plane to either the (110) or (211) MRSS plane at 5 K increases with decreasing strain rate.Extrapolating these observations, it may be expected that the effective MD trajectory should be parallel to the (110) slip plane at very low strain rates, as observed experimentally by most authors [83].
Interestingly, the enthalpy of kink-pair formation, determined with the MTP from the temperature dependence of the CRSS, is very similar to the value obtained with the EAM potential [2].Both values (∼0.03 eV) are much lower than the experimental value, reported by Seeger and Holzwart [80].The discrepancy is thus not necessarily related to limitations of the used interatomic potential.
The developed MTP enables the study of other large-scale defects like 2D dislocation loops as well as vacancy-dislocation interactions and the determination of the Gibbs energy of kink pairs in bcc Nb using hyperdynamics simulations [3] with DFT accuracy.Beyond Nb, the proposed fitting strategy can be applied to obtain accurate MLIPs to model screw dislocations for other bcc elements and alloys.A particularly interesting direction is the application to bcc multi-principal element alloys, such as high entropy alloys, for which the chemical configurational complexity needs to be added to the training set.

Figure 1 .
Figure 1.Typical screw dislocation configurations in the training subset ND4, identified with OVITO [93]: (a) straight screw dislocation; (b) screw dislocation with a single kinkpair and (c) with multiple kink pairs.The X axis is parallel to the [11 2] direction, the Y axis is parallel to the dislocation line ([111] direction) and Z is perpendicular to the (110) slip plane.Red and blue represent screw and edge components, respectively.

Figure 2 .
Figure 2. (a) Energies of atomic configurations in the training set, computed with the MTP, versus the DFT energies (blue dots).The energy root-mean square error (RMSE) is given at the top of the plot.(b) Forces of the atoms in the atomic configurations in the training set, computed with the MTP, versus the DFT forces (blue dots).The forcecomponent RMSE is given at the top of the plot.The black dashed lines represent a perfect fit.

Figure 3 .
Figure 3. DFT energies per atom for the bcc (black empty symbols) and the fcc phases (red empty symbols) as a function of the molar volume.The full lines are the corresponding MTP energies.

Figure 4 .
Figure 4. Energies per atom with respect to the equilibrium bcc energy along the (a) tetragonal, (b) orthorhombic and (c) trigonal deformation paths: empty symbols-DFT calculations, black lines-MTP results.

Figure 6 .
Figure 6.Phonon dispersion relations at 300 K: full circles-experimental data for the longitudinal (L) modes, open circles and full squares-experimental data for the transverse (T) modes; full lines-MTP longitudinal modes, dashed and dotted lines-MTP transverse modes.

Figure 7 .
Figure 7. DFT-MTP relative errors for: (a) properties included in the training set, (b) properties not directly included in the training set.The symbols f (H), f (P) and f (N) denote the phonon frequencies at the H, P and N zone boundaries.

Figure 11 .
Figure 11.Dependence of the critical resolved shear stress (CRSS) on the orientation angle χ (full circles), calculated with the MTP.The full line is a non-linear fit using equation (3), the dashed line is the Schmid law τc/ cos(χ).

Figure 12 .
Figure 12.Dependence of the critical waiting time for the nucleation of the first kink pair on the strain rate at 5 K: (•) (211) PAD models, (o) (110) PAD models.The dashed lines are linear fits through the data points.

Figure 14 .
Figure 14.Initial screw dislocation trajectories for the (211) PAD model at different temperatures.

Figure 15 .
Figure 15.Nucleation and growth of a double kink at 5 K and a strain rate of 5.36 × 10 8 s −1 , shown as a projection perpendicular to the X axis: (a) A just nucleated kink pair with a critical length of ∼2 Å at a time step 69.1 ps and shear stress of 1.57 GPa; (b)-(e) kink pair growth from time step 69.2 to 69.8 ps; (f) the screw dislocation is completely at the next position at a time step of 69.9 ps.The bcc atoms are shown in blue, the core atoms in white.The screw and edge components of the dislocation line are shown in red and blue, respectively.The pink arrows are the relative atomic displacements with respect to the first frame at 68.9 ps.

Figure 16 .
Figure 16.Temperature dependence of the critical resolved shear stress (CRSS): (•) and (■) symbols-SD and PAD results using EAM potential; (o) symbols-PAD results using MTP.The dashed lines are fits with equation (6) in the text.

Figure 17 .
Figure 17.Computational speed of the developed MTP in comparison to several other interatomic potentials for bcc Nb.The speed units are ns of computational time per 1 day of real (physical) time.

Table 1 .
Summary of the training set.All dislocations refer to screw dislocations.Quantities fitted: energy (E), forces (F), stresses (S).

Table 2 .
Basic static properties of bcc Nb obtained by our MTP, DFT calculations, previous MLIPs and experimental measurements.MTP a MTP b SNAP c GAP d DFT * TB e DFT f Exp.

Table 4 .
Comparison of several thermodynamic properties of bcc Nb.

Table 5 .
Static properties of screw dislocations in bcc Nb, predicted with the MTP.

Table 8 .
Orientation of the supercells at different angles χ.