Modelling and Simulation in Materials Science and Engineering Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 (33pp) https://doi.org/10.1088/1361-651X/ad2d68 Moment tensor potential for static and dynamic investigations of screw dislocations in bcc Nb Nikolay Zotov∗, Konstantin Gubaev, Julian Wörner and Blazej Grabowski Institute for Materials Science, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany E-mail: zotov@imw.uni-stuttgart.de Received 30 October 2023; revised 12 February 2024 Accepted for publication 27 February 2024 Published 8 March 2024 Abstract 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 temperat- ure 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 molecu- lar dynamics simulations with near DFT accuracy of properties such as for example the Peierls stress, the critical waiting time for the onset of screw dis- location 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 atommethod potentials thus becomes possible. Keywords: bcc Nb, machine learning, moment tensor potential, molecular dynamics, screw dislocation mobility, kink pairs ∗ Author to whom any correspondence should be addressed. Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. © 2024 The Author(s). Published by IOP Publishing Ltd 1 https://doi.org/10.1088/1361-651X/ad2d68 https://orcid.org/0000-0002-6098-4086 https://orcid.org/0000-0003-2612-8515 https://orcid.org/0000-0003-4281-5665 mailto:zotov@imw.uni-stuttgart.de http://crossmark.crossref.org/dialog/?doi=10.1088/1361-651X/ad2d68&domain=pdf&date_stamp=2024-3-8 https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 1. 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 kink- pair formation were recently performed for the case of bcc Nb with standard and acceler- ated molecular dynamics (MD) utilizing an embedded-atom method (EAM) potential [2–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 bar- rier 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 accur- acy 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 prom- ising new approach in material science for the construction of interatomic potentials, com- bining 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–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–8]. Despite the great recent success, the development of general-purposeMLIPs 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–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 aug- mented with an analytical repulsive potential for radiation damage simulations involving high- energy collisions [12]. A moment tensor potential (MTP) [13] was developed for the predica- tion 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 accur- acy 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 2 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 dis- location 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 c44 elastic constant, the lowest elastic anisotropy (Zener ratio), defined as A= 2c44/(c11 − c12), and the lowest G111 shear modulus (G111 = (c11 − c12 + c44)/3) among all bcc transition metals. Moreover, c44 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 longitud- inal 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 addi- tion, MTP generally shows an excellent trade-off between execution speed and accuracy, com- pared to other MLIPs [25]. 2. Computational methods 2.1. 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–29]. The theoretical description of theMTP 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 amongMTPs even if the potentials exhibit similar root-mean- square (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 config- urations, 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. 3 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 2.2. 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–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 prop- erties is still an open question in theMLIP 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 anMLIP, 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–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 prop- erties. 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 priorimay 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 sym- metry or could not be fitted reliably. The details of each group are described below. 2.2.1. Group A. The training subsets in group A aim at an accurate reproduction of the gen- eral 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 4 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Table 1. Summary of the training set. All dislocations refer to screw dislocations. Quantities fitted: energy (E), forces (F), stresses (S). Group Training subset Number of atoms Total cfgs. What to fit Description A NA1 54 7 ES bcc phase, 7 volumes from 9.8 to 27.4 Å3 NA2 32 7 ES fcc phase, 7 volumes from 16 to 22.8 Å3 NA3 16 18 ES Bain deformation path NA4 4 15 ES Orthorhombic deformation path NA5 96 21 ES Trigonal deformation path NA6 54 113 EFS DFT snapshots at 300 K B NB1 16 5 ES c44 volume-preserving deformation NB2 16 5 ES C’ volume-preserving deformation NB3 16 9 ES c11 deformation mode NB4 16 9 ES c44 deformation mode NB5 16 9 ES c11+c44 deformation mode C NC1 64 21 EFS (110) stacking-fault energy curve along the [111] direction NC2 72 21 EFS (211) stacking-fault energy curve along the [111] direction D ND1 127 15 EFS Configurations with single vacancy ND2 336–480 15 EFS Configurations with (100), (110), (111) or (211) free surfaces ND3 384–600 6 EFS Configurations with two dislocations with (110) or (211) slip planes ND4 216 11 ES Single-dislocation configurations; straight dislocations and dislocations with single and multiple kink-pairs ND5 231 9 ES Dipole models in quadrupole setting along the (110) NEB profile paths). The paths connect the bcc with the fcc, sc and bct structures through a series of non- equilibrium 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 (∆Efcc−bcc) is crucial, because recently it has been shown that the screw disloca- tion core structure is governed by ∆Efcc−bcc ([34] and references therein). The orthorhombic 5 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 shear- ing 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 homo- geneous 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 run- ning 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. 2.2.2. Group B. The stress and dilatation fields of 1 2 〈111〉 screw dislocations depend sens- itively on the elastic constants [36]. The nucleation energy of a kink-pair is proportional to G111b2/2π [2, 37], where b is the length of the Burgers vector. Correspondingly, the training subsets in group B target the correct prediction of the elastic constants. Training subsets NB1 and NB2 target the reproduction of the c44 and C ′ = (c11 − c12)/2 elastic constants, which measure the resistance against different types of shear under constant volume, while NB3 determines the c11 constant. Training subsets NB4 and NB5 are included in order to better con- strain c11 and c44 at larger deformations. The corresponding deformation modes are described in [38]. 2.2.3. Group C. 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 previ- ously for fitting the GAP [12] and the MTP [13] for bcc Nb. 2.2.4. 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 prop- erly 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 6 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 1. Typical screw dislocation configurations in the training subset ND4, identified with OVITO [93]: (a) straight screw dislocation; (b) screw dislocationwith a single kink- pair and (c) with multiple kink pairs. The X axis is parallel to the [112̄] 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. in ND5 represent dipoles with 1b length of the dislocation line in a special quadrupole arrange- ment (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–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 bar- rier, 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–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 screw- dislocation configurations with and without kink pairs. 2.3. 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 Nbwith the valence con- figuration of 4p64d45 s1 (11 valence electrons), distributed with VASP. Except for the energy of the isolated Nb atom, we used non-spin-polarized calculations. We employed the gener- alized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) exchange- correlation functional. On the basis of convergence tests for the cohesive energy, the plane- wave 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 7 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 integ- ration 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. 2.4. 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 pro- gram 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 gradi- ent 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 || [112̄], Y || [111] and Z || [11̄0]. For the dislocation models with (211) MRSS plane, the orientation of the supercell was X || [011̄], Y || [111] and Z || [21̄1̄] for the twinning setup and X || [11̄0], Y || [111] and Z || [112̄] for the anti-twinning setup. For the dislocation models with (123) MRSS plane, the orientation of the supercell was X || [5̄41], Y || [111] and Z || [123̄]. Periodic boundary conditions along Y and free surfaces along X and Z were used in the case of single-dislocation configurations. For the calculation of the elastic constants we used the energy-strain method [38] and strain deformations from −5% to +5% with a step size of 1%. For the static calculations of the SFE curves along the [111] direction, we used periodic boundary conditions in all three directions and the following orientations: X || [001], Y || [110] and Z || [1̄10] for the (110) γ-curves; X || [11̄0], Y || [111̄] and Z || [112] for the (211) γ-curves. For the static calculations of the Peierls barrier between neighbouring equilib- rium easy-core positions we used a skewed supercell containing 231 atoms and a pair of screw dislocations with opposite Burgers vectors (dipole) in a quadrupole arrangement [41]. The basis vectors of the supercell were a1 = 7e1, a2 = 1 2 (7e1 + 11e2 + 1e3), a3 = 1e3, where e1, e2 and e3 are the basis vectors along the [112̄], [1̄10] and the [111] directions, respectively. The same supercell was also used for some of the Peierls stress calculations. The single-dislocation configurations in training subset ND4 containing 216 atoms were oriented as X || [112̄], Y || [111] and Z || [11̄0]. For the static calculations of the Peierls stress using large supercells, energy and force minimization was first performed in the unstressed state, using the conjug- ate 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 (Evac). 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×106 8 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 force- component RMSE is given at the top of the plot. The black dashed lines represent a perfect fit. time steps for equilibration and 3×106 times steps for the determination of the phonon fre- quencies. 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 disper- sion curves were calculated from the dynamical matrix using the program PHANA [56]. The mean square displacement 〈u2〉 at different temperatures was calculated usingMD simulations (1×104 time steps for equilibration and 5×104 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: α= (1/3B)dp/dT, (1) where B is the bulk modulus and p is the pressure. MD simulations using a 9× 9× 9 con- ventional 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). 2.5. Validation of the MTP An accurate and unbiased agreement between the DFT andMTP energies and forces is demon- strated 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 9 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Table 2. Basic static properties of bcc Nb obtained by our MTP, DFT calculations, pre- vious MLIPs and experimental measurements. Constant MTP∗ MTPa MTPb SNAPc GAPd DFT∗ TBe DFTf Exp. Lattice Parameter ao(A) 3.322(1) 3.323 3.316 3.327 3.308 3.322 3.314 3.324 3.303g Ecoh (eV/atom) 6.37 7.004 6.92 6.91 7.523h Elastic Constants c11 (GPa) 247.0(5) 269 222.5(3)k 266 243 244.4(5) 249.0 253i c12 (GPa) 139.0(1.0) 127 153(1)k 142 137 134.8(1.6) 135.4 134i c44 (GPa) 18.5(1.5) 16 25.8(3)k 20 13 17.0(5a) 18.1 30.9i Zener Ratio A 0.34 0.22 0.74 0.32 0.24 0.31 0.32 0.52 Bulk Modulus B (GPa) 175.0(7) 174 176 183 172 171 152 173 174 Shear Moduli G111 (GPa) 42.1(6) 52.6 31.7 47.3 39.6 42.2 43.9 49.9 G ′ l 4 (GPa) 32.9(6) 33.1 31.1 36.5 26.2 31.5 33.2 45.5 Single Vacancy Evac unrelaxed (eV) 2.82(1) 1.78k 3.01 3.53 2.6–3.1j Evac relaxed (eV) 2.51(2) 1.65k 2.85 2.78 3.23 ∗Present results. a [14]. b [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: G ′ 4 = 3c44/(1+ 2A). 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 Ecoh of the bcc and fcc phases (see figure 3) as well as the energy difference between the fcc and the bcc structures∆Efcc−bcc. The elastic constants are also in good agreement with our and previous [58] DFT calculations. Similar results for c11 and c12 are obtained using GAP [12], but GAP [12] significantly underestimates c44. Both theMTP, developed in [14], and the SNAP from [11] overestimate c11, while the MTP, from [13], underestimates c11 and overestimates c12 as well as the Zener ratio. All MLIPs for bcc Nb underestimate the c44 constant, which 10 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 corres- ponding MTP energies. 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. reflects the inability of the underlying DFT methods to correctly predict c44, 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 accur- ately 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 cal- culations (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 max- imum at a displacement 0.5b, are compared in table 3 with previous results. The unrelaxed 11 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 5. Stacking-fault energies along the [111] direction for (a) the (110) slip plane and (b) the (211) slip plane: open symbols—DFT calculations; black lines—MTP results. Table 3. Unrelaxed and relaxed energies of unstable stacking faults (in eVÅ−2) along the (110), (211) and (123) slip planes. SFE MTP∗a DFT∗a DFT1a SNAP2a MTP∗b MTP3b SNAP4b DFT5b (110) 0.049 0.050 0.057 0.056 0.044 0.040 0.057 0.042 (211) 0.058 0.058 0.064 0.052 0.067 0.048 (123) 0.057 0.051 0.064 0.048 ∗Present results. a Unrelaxed. b Relaxed. 1 [63]. 2 [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 sim- ulations 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 vibra- tional properties of bcc Nb, not represented in the training set, are compared with experiments in table 4. The Grüneisen parameter is calculated as αVmolB/Cp , where Vmol is the molar volume and Cp 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 12 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 trans- verse (T) modes; full lines—MTP longitudinal modes, dashed and dotted lines—MTP transverse modes. Table 4. Comparison of several thermodynamic properties of bcc Nb. Property MTP∗ GAP1 TB2 Exp Linear expansion coefficient α at 300 K (×10−6 K−1) 7.6(6) 8.5 5.3 7.1–7.2a Grüneisen parameter at 300 K 1.75(5) 2.0 1.6b Mean-square displacement at 300 K (Å2) 0.0066(1) 0.0039 0.0068(11)c at 90 K (Å2) 0.0019(1) 0.0018d ∗Present results. 1 [12]. 2 TB (tight binding) [66]. a [67]. b [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. c44 versus c12). 3. 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 13 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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. thoroughly validatedMTP (section 2.5). Sections 3.1 and 3.2 contain static and dynamic screw- dislocation properties, respectively, while section 3.3 presents CPU timings for the large-scale dislocation simulations with the MTP. 3.1. 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–44, 72] and empirical potentials [2, 40, 72, 73] for bcc metals. The energy of the easy core Ecore was determined using (110) dipole models of increasing size along X and the equation [72]: ∆E= Ecore +Ks ( b2/4π ) ln(d/rcore) , (2) 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, Ks is a shear modulus, depending on the modified elastic compliances [2], d is the distance between the two dislocations along the X axis, rcore 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 (rcore ∼ 1.633b) was used in equation (2). A linear dependence of ∆E versus ln(d/rcore) 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 Ks = 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 14 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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. Table 5. Static properties of screw dislocations in bcc Nb, predicted with the MTP. Property MTP∗ DFT∗ DFT1 DFT2 Dislocation core NDa NDa NDa Easy core energy, eV/b 0.22(12) 0.201 Peierls barrier 31.5 32.7 29.6 36.0 (110) slip plane, meV/b Peierls barrier profile SHb SHb SHb SHb ∗Present results. 1 [41]. 2 [42]. a ND: non-degenerate core. b SH: single-hump profile. 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. 15 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 9. Peierls barrier: filled symbols—DFT calculations, full line—MTP results. Table 6. Comparison of Peierls stresses (in GPa), calculated with our MTP, with previ- ously published DFT and MLIP results. Potential Model Slip plane Sense of shearing Number of atoms MTP∗ SNAP1 DFT∗ DFT2 Dipole (110) + 231 1.38 1.04 0.74 PAD (110) + 49 980 1.42 0.84 — 1.42 0.83 PAD (211) + 49 980 1.47 1.45 — 2.48 2.10 PAD (123) + 51 072 1.48 1.18 — 2.16 0.54 ∗Present results. 1 [64]. 2 [42]. 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 res- ults 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 a2 basis vector, a ′ 2 = a2 + xe3 (0< x< 1) [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 xcrit (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 xcrit. 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 16 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 10. Dependence of the (211) Peierls stresses for large atomic configurations, con- taining 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 squares— anti-twining orientation with negative sense of shearing (AT−), empty circles—twining orientation with negative sense of shearing (T−), empty triangles—anti-twinning ori- entation 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 c11 and c12 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 dislo- cation 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–44]) to investigate the effect of supercell size along the glide direction (Lx) 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 condi- tions. 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 Lx " 100 Å, regardless of the boundary conditions used. For the (211) slip plane and the SD models, the Peierls stresses are well converged for Lx " 280 Å (see figure 10). Figure 10 also demonstrates that the Peierls stresses, calculated with our MTP, comply very well with the important sym- metry requirement [40] that the stresses for the twinning configuration and a positive sense of 17 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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(χ). 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 Lx " 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 (101̄) slip plane. The value χ =−30◦ corres- ponds to the (2̄11) twinning slip plane and χ = 30◦ to the anti-twinning (1̄1̄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 correspond- ing 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 (1̄10) 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] CRSS(χ) = [τc − τ∗ (a2 sin(2χ)+ a3 cos(2χ+π/6))]/[cos(χ)+ a1 cos(χ +π/3)], (3) 18 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al where τ∗ is the non-glide stress acting perpendicular to the screw dislocation, τc is the CRSS at χ = 0◦ and a1, a2 and a3 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, cal- culated 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 simultan- eously a non-degenerate dislocation core and an accurate Peierls potential, which is not the case for many other interatomic potentials includingMLIPs 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 observa- tions 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. 3.2. 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 constant- strain 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. 3.2.1. Strain rate effects. The extensive literature on the mechanical properties of bcc Nb, studied experimentally at different strain rates up to 1.1×103 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 19 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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. meaningful shear strains in reasonable computational time. The developed MTP was tested and shows good robustness for strain rates up to 5×109 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 practic- ally constant with increasing strain rate suggesting that different mechanisms of screw dislo- cation mobility might be operational at low and high strain rates. The CRSS τc is commonly expressed as τc = Gγc, (4) 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 tc for the onset of screw dislocation movement should decrease inversely with increasing strain rate at constant CRSS (taking into account that γc = γ̇tc) tc = t0/γ̇, (5) where t0 is a fit parameter and γ̇ is the strain rate. Indeed, plotting tc versus the strain rate on a log-log scale, a linear increase of tc 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. 3.2.2. 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 20 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Table 7. Characteristics of the screw dislocation trajectories at 5 K. MRSS plane Pass∗ Elementary slip planes Effective trajectory (11̄0) 0 (12̄1)+(1̄1̄2) Initially parallel to (01̄1), followed by multiple cross-slip events "1 Not well defined Wavy trajectory parallel to (11̄0) (21̄1̄) 0 (21̄1̄)+(12̄1) Initially parallel to (1̄10), followed by wavy trajectory parallel to (21̄1̄) "1 (21̄1̄)+(12̄1)+(1̄1̄2) Parallel to (21̄1̄) ∗Pass 0: The screw dislocation reaches the right side of the supercell along the X axis. Pass n: The screw dislocation transverses n times (n " 1) the supercell along the X axis. different models of kink-pair nucleation, and concluded that their data is in quantitative agree- ment 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, inde- pendent 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 (01̄1) 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 (11̄0) 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 (01̄1) slip plane to an effective glide along the (11̄0) 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 (12̄1) 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 (1̄10) 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 (21̄1̄) 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 (1̄10) slip plane to an effective glide along the (21̄1̄) 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 (1̄10) 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 (21̄1̄) 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 predomin- antly parallel to alternating {211} symmetry-equivalent planes: (12̄1)+(1̄1̄2) in the case of the (11̄0) models (figure 13(a)) and (21̄1̄)+(12̄1) in the case of the (21̄1̄) models (figure 13(b)). In 21 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 13. (a) Initial and (b) subsequent screw trajectory for the (110) PAD model; (c) Initial and (d) subsequent screw dislocation trajectory for the (211) PAD model. The strain rate is 6.18×108 s−1 for the (211) PAD model and 5.36×108 s−1 for the (110) PAD model. the subsequent traverses of the supercell the elementary jumps are approximately parallel to (11̄0) in the case of the (11̄0) models (figure 13(c)) and (21̄1̄)+(12̄1)+(1̄1̄2) in the case of the (21̄1̄) models (figure 13(d)). 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 dir- ection, 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]. 22 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 14. Initial screw dislocation trajectories for the (211) PAD model at different temperatures. 3.2.3. Temperature dependence of the CRSS. Already at 5 K the onset of the screw dislo- cation 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–88]) obtained with EAM potentials. To assess the thermally-activated regime of kink-pair nucleation and growth with the developed MTP over a broader temper- ature range, the CRSS was determined as a function of temperature, utilizing a strain rate of 6.18×108 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 con- stant athermal stress level τ a of about 0.55 GPa for T" 400 K. The presence of a ‘hump’ in 23 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Figure 15. Nucleation and growth of a double kink at 5 K and a strain rate of 5.36×108 s−1, shown as a projection perpendicular to the X axis: (a) A just nucle- ated 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 dislo- cation 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 disloca- tion 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. 24 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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. the CRSS-temperature relationship is predicted by the theory of Seeger [23]. In the thermally- activated regime, the temperature dependence of the CRSS can be phenomenologically written in the form [2, 89] τc = τa+ τ0 { 1− [ − kBT ∆Hkp (0) ln ( γ̇ γ̇0 )] 1 q { 1 p , (6) where ∆Hkp(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× 109 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, ∆Hkp(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, ∆Hkp(0) = 0.03± 0.01 eV, p= 0.5 and q= 1.45 [2]. The differences in 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). 25 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 neces- sarily 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 (11̄0)+ (01̄1) slip planes in the case of (211) PAD models and along (01̄1)+ (1̄01) 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,∆Hkp, 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 real- istic 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 bar- rier, using NEB calculations at 0 K. These parameters were then used to calculate numeric- ally the enthalpy of kink pair formation using a linearized version of the line-tension func- tional. Specifically, for bcc Nb they reported a value (1.28 eV), which is almost twice lar- ger than the experimental value (0.64 eV [80]). We expect the here presented dynamic cal- culations at the relevant temperatures with a verified MLIP to better reflect the experimental conditions. The remaining discrepancy between experiment and simulations could originate in the dif- ferent 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 com- bined 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. 3.3. Computational speed We compare here the computational performance of the developed MTP with classical poten- tials and other MLIPs. Specifically, we compare with the EAM potential used in our previous MD simulations [2–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. 26 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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. 4. 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 mechan- ical, 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 train- ing 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 repro- duces 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 (1̄01) 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 ourMTP (level 16) is relatively simple compared to other MLmodels, 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 27 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al 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 prop- erties. 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 temper- atures. 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] (unlikeMTP) render theMTPmore 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×107 to 2.5×109 s−1, leading to an increase of the critical waiting time for the onset of the screw dislocation move- ment 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 find- ing 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 tem- perature 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 developedMTP 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 configur- ational complexity needs to be added to the training set. 28 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al Data availability statement The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors. Acknowledgment This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 865855). Support by the Stuttgart Center for Simulation Science (SimTech) is gratefully acknowledged. The authors also acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 40/575-1 FUGG (JUSTUS 2 cluster). The authors thankMHodapp and A Shapeev for provid- ing their MTP for bcc Nb. N Z thanks YMishin for providing a model of bcc Ta in quadrupole setting. Conflict of interest No potential conflict of interest was reported by the authors. Appendix The orientations of the supercells used for the analysis of the non-Schmid behavior are given in table 8 below. For brevity, only the orientations for χ > 0 are given. The orientations for χ < 0 are obtained by exchanging the h and lMiller indices. Table 8. Orientation of the supercells at different angles χ. Angle χ X axis Z axis 30 (1 1̄ 0) (1 1 2̄) 26 (1 1̄4 13) (9̄ 4 5) 23 (1 8̄ 7) (5̄ 2 3) 19 (1 5̄ 4) (3̄ 1 2) 14 (2 7̄ 5) (4̄ 1 3) 9 (4 1̄1 7) (6̄ 1 5) 7 (6 1̄5 9) (8̄ 1 7) 5 (8 1̄9 11) (1̄0 1 9) ORCID iDs Nikolay Zotov https://orcid.org/0000-0002-6098-4086 Konstantin Gubaev https://orcid.org/0000-0003-2612-8515 Blazej Grabowski https://orcid.org/0000-0003-4281-5665 29 https://orcid.org/0000-0002-6098-4086 https://orcid.org/0000-0002-6098-4086 https://orcid.org/0000-0003-2612-8515 https://orcid.org/0000-0003-2612-8515 https://orcid.org/0000-0003-4281-5665 https://orcid.org/0000-0003-4281-5665 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al References [1] Hirth J P and Lothe J 2015 Theory of Dislocations (Wiley) [2] Zotov N and Grabowski B 2021 Molecular dynamics simulations of screw dislocation mobility in bcc Nb Model. Simul. Mater. Sci. Eng. 29 085007 [3] Grabowski B and Zotov N 2021 Thermally-activated dislocation mobility in bcc metals: an accel- erated molecular dynamics study Comput. Mater. Sci. 200 110804 [4] Zotov N and Grabowski B 2022 Entropy of kink pair formation on screw dislocations: an acceler- ated molecular dynamics study Model. Simul. Mater. Sci. Eng. 30 065004 [5] Farkas D and Jones C 1996 Interatomic potentials for ternary Nb - Ti - Al alloys Model. Simul. Mater. Sci. Eng. 4 23–32 [6] Behler J 2016 Perspective: machine learning potentials for atomistic simulations J. Chem. Phys. 145 170901 [7] Mueller T, Hernandez A and Wang C 2020 Machine learning for interatomic potential models J. Chem. Phys. 152 050902 [8] Mishin Y 2021 Machine-learning interatomic potentials for materials science Acta Mater. 214 116980 [9] Behler J and Csányi G 2021 Machine learning potentials for extended systems: a perspective Eur. J. B 94 142–53 [10] Miksch A M, Morawietz T, Kästner J, Urban A and Artrith N 2021 Strategies for the construction of machine-learning potentials for accurate and efficient atomic-scale simulationsMach. Learn. Sci. Technol. 2 031001 [11] Li X G, Chen C, Zheng H, Zuo Y and Ong S 2020 Complex strengthening mechanisms in the NbMoTaW multi-principal element alloy npj Comput. Mater. 6 70 [12] Byggmästar J, Nordlund K and Djurabekova F 2020 Gaussian approximation potentials for body- centered-cubic transition metals Phys. Rev. Mater. 4 093802 [13] Hodapp M and Shapeev A 2021 Machine-learning potentials enable predictive and tractable high- throughput screening of random alloys Phys. Rev. Mater. 5 113802 [14] KwonH, ShigaM, Kimizuka H and Oda T 2023 Accurate description of hydrogen diffusivity in bcc metals using machine-learning moment tensor potentials and path-integral methods Acta Mater. 247 118739 [15] Jung J H, Srinivasan P, Forslund A and Grabowski B 2023 High-accuracy themodynamic properties to themelting point from ab initio calculations aided bymachine-learning potentials npj Comput. Mater. 9 1–12 [16] Forslund A, Jung J H, Srinivasan P and Grabowski B 2023 Thermodynamic properties on the homo- logous temperature scale from direct upsampling: understanding electron-vibration coupling and thermal vacancies in bcc refractory metals Phys. Rev. B 107 174309 [17] Srinivasan P, Demuriya D, Grabowski B and Shapeev A 2024 Electronic moment tensor potentials include both electronic and vibrational degrees of freedom npj Comput. Mater. 10 41 [18] Nagasako N, Jahnatek M, Asahi R and Hafner J 2010 Anomalies in the response of V, Nb and Ta to tensile and shear loading: ab initio density functional calculations Phys. Rev. B 81 094108 [19] Ashkenazi J, Dacorongna M, Peter M, Talmor Y and Walker E 1978 Elastic constants in Nb-Zr alloys from zero temperature to the melting point: experiment and theory Phys. Rev. B 18 4120– 31 [20] Nakagawa Y and Woods A D B 1963 Lattice dynamics of niobium Phys. Rev. Lett. 11 271–4 [21] Sharp R I 1969 The lattice dynamics of niobium I. Measurements of the phonon frequencies J. Phys. C: Solid State Phys 2 421–31 [22] Woods A D 1964 Lattice dynamics of tantalum Phys. Rev. A 136 781–3 [23] Seeger A 2002 Peierls barriers, kinks and flow stress: recent progress Z. Met. 93 760–77 [24] Shapeev A V 2016 Moment tensor potentials: a class of systematically improvable interatomic potentials Multiscale Model. Simul. 14 1153–73 [25] Zuo Y et al 2020 Performance and cost assessment of machine learning interatomic potentials J. Phys. Chem. A 124 731–45 [26] Gubaev K, Podryabinkin E V, Hart G L and Shapeev A V 2019 Accelerating high-throughput searches for new alloys with active learning of interatomic potentials Comput. Mater. Sci. 156 148–56 30 https://doi.org/10.1088/1361-651X/ac2b02 https://doi.org/10.1088/1361-651X/ac2b02 https://doi.org/10.1016/j.commatsci.2021.110804 https://doi.org/10.1016/j.commatsci.2021.110804 https://doi.org/10.1088/1361-651X/ac7ac9 https://doi.org/10.1088/1361-651X/ac7ac9 https://doi.org/10.1088/0965-0393/4/1/004 https://doi.org/10.1088/0965-0393/4/1/004 https://doi.org/10.1063/1.4966192 https://doi.org/10.1063/1.4966192 https://doi.org/10.1063/1.5126336 https://doi.org/10.1063/1.5126336 https://doi.org/10.1016/j.actamat.2021.116980 https://doi.org/10.1016/j.actamat.2021.116980 https://doi.org/10.1140/epjb/s10051-021-00156-1 https://doi.org/10.1140/epjb/s10051-021-00156-1 https://doi.org/10.1088/2632-2153/abfd96 https://doi.org/10.1088/2632-2153/abfd96 https://doi.org/10.1038/s41524-020-0339-0 https://doi.org/10.1038/s41524-020-0339-0 https://doi.org/10.1103/PhysRevMaterials.4.093802 https://doi.org/10.1103/PhysRevMaterials.4.093802 https://doi.org/10.1103/PhysRevMaterials.5.113802 https://doi.org/10.1103/PhysRevMaterials.5.113802 https://doi.org/10.1016/j.actamat.2023.118739 https://doi.org/10.1016/j.actamat.2023.118739 https://doi.org/10.1038/s41524-022-00956-8 https://doi.org/10.1038/s41524-022-00956-8 https://doi.org/10.1103/PhysRevB.107.174309 https://doi.org/10.1103/PhysRevB.107.174309 https://doi.org/10.1038/s41524-024-01222-9 https://doi.org/10.1038/s41524-024-01222-9 https://doi.org/10.1103/PhysRevB.81.094108 https://doi.org/10.1103/PhysRevB.81.094108 https://doi.org/10.1103/PhysRevB.18.4120 https://doi.org/10.1103/PhysRevB.18.4120 https://doi.org/10.1103/PhysRevB.18.4120 https://doi.org/10.1103/PhysRevLett.11.271 https://doi.org/10.1103/PhysRevLett.11.271 https://doi.org/10.1088/0022-3719/2/3/306 https://doi.org/10.1088/0022-3719/2/3/306 https://doi.org/10.1103/PhysRev.136.A781 https://doi.org/10.1103/PhysRev.136.A781 https://doi.org/10.3139/146.020760 https://doi.org/10.3139/146.020760 https://doi.org/10.1137/15M1054183 https://doi.org/10.1137/15M1054183 https://doi.org/10.1021/acs.jpca.9b08723 https://doi.org/10.1021/acs.jpca.9b08723 https://doi.org/10.1016/j.commatsci.2018.09.031 https://doi.org/10.1016/j.commatsci.2018.09.031 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al [27] Gubaev K, Zaverkin V, Srinivasan P, Duff A I, Kästner J and Grabowski B 2023 Performance of two complementary machine-learned potentials in modelling chemically complex systems npj Comput. Mater. 9 129 [28] Rosenbrock C W, Gubaev K, Shapeev A V, Pártay L B, Bernstein N, Csányi G and Hart G L 2021 Machine-learned interatomic potentials for alloys and alloy phase diagrams npj Comput. Mater. 7 24 [29] Shapeev A V, Podryabinkin E V, Gubaev K, Tasnádi F and Abrikosov I A 2020 Elinvar effect in β-Ti simulated by on-the-fly trained moment tensor potential New J. Phys. 22 113005 [30] Novikov I S, Gubaev K, Podryabinkin E V and Shapeev A V 2020 The MLIP package: moment tensor potentials with MPI and active learning Mach. Learn. Sci. Technol. 2 025002 [31] Bianchini F, Glielmo A, Kermode J R and De Vita A 2019 Enabling qm-accurate simulation of dis- location motion in γ-Ni andα-Fe using a hybrid multiscale approach Phys. Rev. Mater. 3 043605 [32] Podryabinkin E V, Tikhonov E V, Shapeev A V and Oganov A R 2019 Accelerating crystal struc- ture prediction by machine-learning interatomic potentials with active learning Phys. Rev. B 99 064114 [33] Paidar V, Wang L G, Sob M and Vitek V 1999 A study of the applicability of many-body central force potentials in NiAl and TiAlModel. Simul. Mater. Sci. Eng. 7 369–81 [34] Wang R, Ma X, Zhang L, Wang H, Srolovitz D J, Wen T and Wu Z 2022 Classical and machine learning interatomic potentials for bcc vanadium Phys. Rev. Mater. 6 113603 [35] LuoW, RoundyD, CohenML andMorris JW 2002 Ideal strength of bccmolybdenum and niobium Phys. Rev. B 66 094110 [36] Chou Y and Mitchell T 1967 Stress and dilatation fields of the 〈111〉 dislocation in cubic crystals J. Appl. Phys. 38 1535–40 [37] Koizumi H, Kirchner H and Suzuki T 1993Kink pair nucleation and critical shear stressActaMetal. Mater. 41 3483–93 [38] Łopuszyński M and Majewski J A 2007 Ab initio calculations of third-order elastic constants and related properties for selected semiconductors Phys. Rev. B 76 045202 [39] Vítek V 1968 Intrinsic stacking faults in body-centred cubic crystals Phil. Mag. 18 773–86 [40] Duesbery M and Vitek V 1998 Plastic anisotropy in b.c.c. transition metals Acta Mater. 46 1481–92 [41] Li J, Wang C Z, Chang J P, Cai W, Bulatov V V, Ho K M and Yip S 2004 Core energy and Peierls stress of a screw dislocation in bcc molybdenum: a periodic-cell tight-binding study Phys. Rev. B 70 104113 [42] Weinberger C R, Tucker G J and Foiles S M 2013 Peierls potential of screw dislocations in bcc transition metals: predictions from density functional theory Phys. Rev. B 87 054114 [43] Dezerald L, Ventelon L, Clouet E, Denoual C, Rodney D and Willaime F 2014 Ab initio modeling of the two-dimensional energy landscape of screw dislocations in bcc transition metals Phys. Rev. B 89 024104 [44] Dezerald L, Proville L, Ventelon L, Willaime F and Rodney D 2015 First-principles prediction of kink-pair activation enthalpy on screw dislocations in bcc transition metals: V, Nb, Ta, Mo, W and Fe Phys. Rev. B 91 094105 [45] Henkelman G and Jonsson H 2000 Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points J. Chem. Phys. 113 9978–85 [46] Hirel P 2015 Atomsk: a tool for manipulating and converting atomic data files Comput. Phys. Commun. 197 212–9 [47] Goryaeva AM, Dérés J, Lapointe C, Grigorev P, Swinburne T D, Kermode J R, Ventelon L, Baima J and Marinica M-C 2021 Efficient and transerable machine learning potentials for the simulation of crystal defects in bcc Fe and W Phys. Rev. Mater. 5 103803 [48] Szlachta W J, Bartok A P and Csanyi G 2014 Accuracy and transferability of Gaussian approxim- ation potential models for tungsten Phys. Rev. B 90 104108 [49] Lin Y S, Pun G P P and Mishin Y 2022 Development of a physically-informed neural network interatomic potential for tantalum Comput. Mater. Sci. 205 111180 [50] Kresse G and Furthmüller J 1996 Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set Comput. Mater. Sci. 6 15–50 [51] Monkhorst H J and Pack J D 1976 Special points for Brillouin-zone integrations Phys. Rev. B 13 5188–92 [52] Methfessel M and Paxton A T 1989 High-precision sampling for Brillouin-zone integration in metals Phys. Rev. B 40 3616–21 31 https://doi.org/10.1038/s41524-023-01073-w https://doi.org/10.1038/s41524-023-01073-w https://doi.org/10.1038/s41524-020-00477-2 https://doi.org/10.1038/s41524-020-00477-2 https://doi.org/10.1088/1367-2630/abc392 https://doi.org/10.1088/1367-2630/abc392 https://doi.org/10.1088/2632-2153/abc9fe https://doi.org/10.1088/2632-2153/abc9fe https://doi.org/10.1103/PhysRevMaterials.3.043605 https://doi.org/10.1103/PhysRevMaterials.3.043605 https://doi.org/10.1103/PhysRevB.99.064114 https://doi.org/10.1103/PhysRevB.99.064114 https://doi.org/10.1088/0965-0393/7/3/306 https://doi.org/10.1088/0965-0393/7/3/306 https://doi.org/10.1103/PhysRevMaterials.6.113603 https://doi.org/10.1103/PhysRevMaterials.6.113603 https://doi.org/10.1103/PhysRevB.66.094110 https://doi.org/10.1103/PhysRevB.66.094110 https://doi.org/10.1063/1.1709719 https://doi.org/10.1063/1.1709719 https://doi.org/10.1016/0956-7151(93)90228-K https://doi.org/10.1016/0956-7151(93)90228-K https://doi.org/10.1103/PhysRevB.76.045202 https://doi.org/10.1103/PhysRevB.76.045202 https://doi.org/10.1080/14786436808227500 https://doi.org/10.1080/14786436808227500 https://doi.org/10.1016/S1359-6454(97)00367-4 https://doi.org/10.1016/S1359-6454(97)00367-4 https://doi.org/10.1103/PhysRevB.70.104113 https://doi.org/10.1103/PhysRevB.70.104113 https://doi.org/10.1103/PhysRevB.87.054114 https://doi.org/10.1103/PhysRevB.87.054114 https://doi.org/10.1103/PhysRevB.89.024104 https://doi.org/10.1103/PhysRevB.89.024104 https://doi.org/10.1103/PhysRevB.91.094105 https://doi.org/10.1103/PhysRevB.91.094105 https://doi.org/10.1063/1.1323224 https://doi.org/10.1063/1.1323224 https://doi.org/10.1016/j.cpc.2015.07.012 https://doi.org/10.1016/j.cpc.2015.07.012 https://doi.org/10.1103/PhysRevMaterials.5.103803 https://doi.org/10.1103/PhysRevMaterials.5.103803 https://doi.org/10.1103/PhysRevB.90.104108 https://doi.org/10.1103/PhysRevB.90.104108 https://doi.org/10.1016/j.commatsci.2021.111180 https://doi.org/10.1016/j.commatsci.2021.111180 https://doi.org/10.1016/0927-0256(96)00008-0 https://doi.org/10.1016/0927-0256(96)00008-0 https://doi.org/10.1103/PhysRevB.13.5188 https://doi.org/10.1103/PhysRevB.13.5188 https://doi.org/10.1103/PhysRevB.40.3616 https://doi.org/10.1103/PhysRevB.40.3616 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al [53] Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 117 1–19 [54] Kong L T 2009 Implementation of Green’s function molecular dynamics: an extension to lammps Comput. Phys. Commun. 180 1004–10 [55] Kong L T 2011 Phonon dispersion measured directly from molecular dynamics simulations Comput. Phys. Commun. 182 2201–7 [56] Kong L T 2020 Phana (Shanghai Jiao Tong University) (available at: https://github.com/lingtikong/ phana) [57] Söderlind P, Yang L H, Moriarty J A and Wills J M 2000 First-principles formation energies of monovacancies in bcc transition metals Phys. Rev. B 61 2579–86 [58] Xu S, Su Y, Smith L and Beyerlein I 2020 Frank-read source operation in six body-centered cubic refractory metals J. Mech. Phys. Solids 141 104017 [59] Roberge R 1975 Lattice parameter of niobium between 4.2 and 300 K J. Less Common Metal. 40 161–4 [60] Kittel C 1996 Introduction to Solid State Physics 7th edn (Wiley) [61] Carroll K J 1965 Elastic constants of niobium from 4.2◦ to 300◦ K J. Appl. Phys. 36 3689–90 [62] Ullmaier H (ed) 1991 Atomic Defects in Metals (Landolt-BöRnstein - Group III Condensed Matter vol 25 (Springer) [63] Frederiksen S L and Jacobsen K W 2003 Density functional theory studies of screw dislocation core structures in bcc metals Phil. Mag. 83 365–75 [64] Wang X, Xu S, Jian W R, Li X G, Su Y and Beyerlein I 2021 Generalized stacking fault ener- gies and Peierls stresses in refractory body-centered cubic metals from machine learning-based interatomic potentials Comput. Mater. Sci. 192 110364 [65] White G K 1988 The heat capacity of transition metals at high temperatures Physica B+C 149 255–60 [66] Lekka C, Bernstein N, Papaconstantopoulos D and Mehl M 2009 Properties of bcc metals by tight- binding total energy simulations Mater. Sci. Eng. B 163 8–16 [67] Thurnay K 1998 Thermal properties of transition metals Sci. Report 6096 (FZKA) [68] Bashir J, Khan Q H and Butt N B 1987 Debye–Waller coefficient of Nb by the elastic neutron diffraction method Acta Cryst. A 43 795–7 [69] Peng L M, Ren G, Dudarev S L and Whelan M J 1996 Debye–Waller factors and absorptive scat- tering factors of elemental crystals Acta Cryst. A 52 456–70 [70] Vítek V, Perrin R C and Bowen D K 1970 The core structure of 1/2 screw dislocations in b.c.c. crystals Phil. Mag. 21 1049–73 [71] Groger R 2017 Ddplot Ver 4.0 (Brno, Institute of Physics of Materials) [72] Ismail-Beigi S and Arias T A 2000 Ab initio study of screw dislocations inMo and Ta: a new picture of plasticity in bcc transition metals Phys. Rev. Lett. 84 1499–502 [73] Ito K and Vitek V 2001 Atomistic study of non-Schmid effects in the plastic yielding of bcc metals Phil. Mag. A 81 1387–407 [74] Shimizu F, Ogata S, Kimizuka H, Kano T, Li J and Kaburaki H 2007 First-principles calculation on screw dislocation core properties in bcc molybdenum J. Earth Simul. 7 17–21 [75] Zhang X C, Cao S, Zhang L J, Yang R and Hu Q M 2022 Unstable stacking fault energy and Peierls stress for evaluating slip system competition in body-centered cubic metals J. Mater. Res. Technol. 22 3413–22 [76] Gröger R, Racherla V, Bassani J L and Vítek V 2008 Multiscale modeling of plastic deformation of molybdenum and tungsten: II. Yield criterion for single crystals based on atomistic studies of glide of 1/2〈111〉 screw dislocations Acta Mater. 56 5412–25 [77] Kamimura Y, Edagawa K and Takeuchi S 2013 Experimental evaluation of the Peierls stresses in a variety of crystals and their relation to the crystal structure Acta Mater. 61 294–309 [78] Starikov S, Grigorev P and Olsson P 2023 Angular-dependent interatomic potential for large-scale atomistic simulations of W-Mo-Nb ternary alloys Comput. Mater. Sci. 233 112734 [79] Maresca F, Dragoni D, Csányi G, Marzari N and Curtin W A 2018 Screw dislocation structure and mobility in body centered cubic Fe predicted by a Gaussian approximation potential npj Comput. Mater. 4 69–76 [80] Seeger A and Holzwarth U 2006 Slip planes and kink properties of screw dislocations in high-purity niobium Phil. Mag. 86 3861–92 32 https://doi.org/10.1006/jcph.1995.1039 https://doi.org/10.1006/jcph.1995.1039 https://doi.org/10.1016/j.cpc.2008.12.035 https://doi.org/10.1016/j.cpc.2008.12.035 https://doi.org/10.1016/j.cpc.2011.04.019 https://doi.org/10.1016/j.cpc.2011.04.019 https://github.com/lingtikong/phana https://github.com/lingtikong/phana https://doi.org/10.1103/PhysRevB.61.2579 https://doi.org/10.1103/PhysRevB.61.2579 https://doi.org/10.1016/j.jmps.2020.104017 https://doi.org/10.1016/j.jmps.2020.104017 https://doi.org/10.1016/0022-5088(75)90193-9 https://doi.org/10.1016/0022-5088(75)90193-9 https://doi.org/10.1063/1.1703072 https://doi.org/10.1063/1.1703072 https://doi.org/10.1080/0141861021000034568 https://doi.org/10.1080/0141861021000034568 https://doi.org/10.1016/j.commatsci.2021.110364 https://doi.org/10.1016/j.commatsci.2021.110364 https://doi.org/10.1016/0378-4363(88)90251-3 https://doi.org/10.1016/0378-4363(88)90251-3 https://doi.org/10.1016/j.mseb.2009.04.014 https://doi.org/10.1016/j.mseb.2009.04.014 https://doi.org/10.1107/S0108767387098507 https://doi.org/10.1107/S0108767387098507 https://doi.org/10.1107/S010876739600089X https://doi.org/10.1107/S010876739600089X https://doi.org/10.1080/14786437008238490 https://doi.org/10.1080/14786437008238490 https://doi.org/10.1103/PhysRevLett.84.1499 https://doi.org/10.1103/PhysRevLett.84.1499 https://doi.org/10.1080/01418610108214447 https://doi.org/10.1080/01418610108214447 https://doi.org/10.32131/jes.7.17 https://doi.org/10.32131/jes.7.17 https://doi.org/10.1016/j.jmrt.2022.12.162 https://doi.org/10.1016/j.jmrt.2022.12.162 https://doi.org/10.1016/j.actamat.2008.07.037 https://doi.org/10.1016/j.actamat.2008.07.037 https://doi.org/10.1016/j.actamat.2012.09.059 https://doi.org/10.1016/j.actamat.2012.09.059 https://doi.org/10.1016/j.commatsci.2023.112.723 https://doi.org/10.1016/j.commatsci.2023.112.723 https://doi.org/10.1038/s41524-018-0125-4 https://doi.org/10.1038/s41524-018-0125-4 https://doi.org/10.1080/14786430500531769 https://doi.org/10.1080/14786430500531769 Modelling Simul. Mater. Sci. Eng. 32 (2024) 035032 N Zotov et al [81] Jian W R, Zhang M, Xu S and Beyerlein I J 2020 Atomistic simulations of dynamics of an edge dislocation and its interaction with a void in copper: a comparative study Model. Simul. Mater. Sci. Eng. 28 045004 [82] Croteau J F et al 2020 Effect of strain rate on tensile mechanical properties of high-purity niobium single crystals for SRF applications Mater. Sci. Eng. A 797 140258 [83] Weinberger C R, Boyce B L and Battaile C C 2013 Slip planes in bcc transition metals Int. Mater. Rev. 58 296–314 [84] Bowen D K, Christian J W and Taylor G 1967 Deformation properties of niobium single crystals Can. J. Phys. 45 903–38 [85] Duesbery M S and Foxall R A 1969 A detailed study of the deformation of high purity niobium single crystals Phil. Mag. 20 719–51 [86] Ngan A HW and Wen M 2002 Atomistic simulation of energetics of motion of screw dislocations in bcc Fe at finite temperatures Comput. Mater. Sci. 23 139–45 [87] Chaussion J, Frivel M and Rodney D 2006 The glide of screw dislocations in bcc Fe: atomistic static and dynamic simulations Acta Mater. 54 3407–16 [88] Shinzato S, Wakeda M and Ogata S 2019 An atomistically informed kinetic Monte Carlo model predicting solid solution strengthening of body-centered cubic alloys Int. J. Plast. 122 319–37 [89] Nemat-Nasser S and Guo W 2000 Flow stress of commercially pure niobium over a broad range of temperatures and strain rates Mater. Sci. Eng. A 284 202–10 [90] Takeuchi S, Hashimoto T andMaeda K 1982 Plastic defomation of bcc metal single crystals at very low temperatures Trans. Jpn. Inst. Met. 23 60–69 [91] Lee B J, Baskes M I, Kim H and Cho Y K 2001 Second nearest-neighbor modified embedded atom method potentials for bcc transition metals Phys. Rev. B 64 184102 [92] Thompson A P, Swiler L P, Trott C R, Foiles S M and Tucker G J 2015 Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials J. Comput. Phys. 285 316–30 [93] Stukowski A 2010 Visualization and analysis of atomistic data with OVITO-the open visualization tool Mod. Simul. Mater. Sci. Eng. 18 015012 33 https://doi.org/10.1088/1361-651X/ab8358 https://doi.org/10.1088/1361-651X/ab8358 https://doi.org/10.1016/j.msea.2020.140258 https://doi.org/10.1016/j.msea.2020.140258 https://doi.org/10.1179/1743280412Y.0000000015 https://doi.org/10.1179/1743280412Y.0000000015 https://doi.org/10.1139/p67-069 https://doi.org/10.1139/p67-069 https://doi.org/10.1080/14786436908228040 https://doi.org/10.1080/14786436908228040 https://doi.org/10.1016/S0927-0256(01)00224-5 https://doi.org/10.1016/S0927-0256(01)00224-5 https://doi.org/10.1016/j.actamat.2006.03.044 https://doi.org/10.1016/j.actamat.2006.03.044 https://doi.org/10.1016/j.ijplas.2019.03.004 https://doi.org/10.1016/j.ijplas.2019.03.004 https://doi.org/10.1016/S0921-5093(00)00740-1 https://doi.org/10.1016/S0921-5093(00)00740-1 https://doi.org/10.2320/matertrans1960.23.60 https://doi.org/10.2320/matertrans1960.23.60 https://doi.org/10.1103/PhysRevB.64.184102 https://doi.org/10.1103/PhysRevB.64.184102 https://doi.org/10.1016/j.jcp.2014.12.018 https://doi.org/10.1016/j.jcp.2014.12.018 https://doi.org/10.1088/0965-0393/18/1/015012 https://doi.org/10.1088/0965-0393/18/1/015012 Moment tensor potential for static and dynamic investigations of screw dislocations in bcc Nb 1. Introduction 2. Computational methods 2.1. MTP 2.2. Training set 2.2.1. Group A. 2.2.2. Group B. 2.2.3. Group C. 2.2.4. Group D. 2.3. DFT calculations 2.4. Static and MD simulations 2.5. Validation of the MTP 3. Results and discussion 3.1. Static screw-dislocation properties 3.2. Dynamic screw-dislocation properties 3.2.1. Strain rate effects. 3.2.2. Trajectories and slip planes. 3.2.3. Temperature dependence of the CRSS. 3.3. Computational speed 4. Summary and conclusions Appendix References