# Activation entropy of dislocation glide in body-centered cubic metals from atomistic simulations

Arnaud Allera, Thomas D Swinburne, Alexandra M Goryaeva, Baptiste

Bienvenu, Fabienne Ribeiro, Michel Perez, Mihai-Cosmin Marinica, David

Rodney

► To cite this version:

Arnaud Allera, Thomas D Swinburne, Alexandra M Goryaeva, Baptiste Bienvenu, Fabienne Ribeiro, et al.. Activation entropy of dislocation glide in body-centered cubic metals from atomistic simulations. Nature Communications, 2025, 16 (1), pp.8367. 10.1038/s41467-025-62390-w. hal-05285237

HAL Id: hal-05285237

https://hal.science/hal-05285237v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Activation entropy of dislocation glide in body-centered cubic metals from atomistic simulations

Received: 11 September 2024

Accepted: 21 July 2025

Published online: 24 September 2025

<figure><img src="image_3.png" /></figure>

Check for updates

Arnaud Allera Ⓔ<fnref n="1" />,²,³,⁴✉, Thomas D. Swinburne Ⓔ⁵,⁶, Alexandra M. Goryaeva Ⓔ³, Baptiste Bienvenu Ⓔ³,⁷, Fabienne Ribeiro<fnref n="1" />, Michel Perez Ⓔ⁴, Mihai-Cosmin Marinica Ⓔ³ & David Rodney Ⓔ²✉

The activation entropy of dislocation glide, a key process controlling the strength of many metals, is often assumed to be constant or linked to enthalpy through the empirical Meyer-Neldel law-both of which are simplified approximations. In this study, we take a more direct approach by calculating the activation Gibbs energy for kink-pair nucleation on screw dislocations of two body-centered cubic metals, iron and tungsten. To ensure reliability, we develop machine learning interatomic potentials for both metals, carefully trained on dislocation data from density functional theory. Our findings reveal that dislocations undergo harmonic transitions between Peierls valleys, with an activation entropy that remains largely constant, regardless of temperature or applied stress. We use these results to parameterize a thermally-activated model of yield stress, which consistently matches experimental data in both iron and tungsten. Our work challenges recent studies using classical potentials, which report highly varying activation entropies, and suggests that simulations relying on classical potentials-widely used in materials modeling-could be significantly influenced by overestimated entropic effects.

A classical thermally-activated process in metallurgy is the glide of screw dislocations in body-centered cubic (BCC) metals<fnref n="1" />. Despite its importance, this process has so far been primarily treated phenomenologically², assuming a simple expression for the activation entropy and a linear relation between the applied stress tensor and the critical resolved shear stress at OK. In particular, the interplay between applied stress and temperature, i.e. the stress- and temperature-dependence of the Gibbs activation energy for dislocation glide, remains a theoretical obstacle where assumptions are commonly used without appropriate justification

Dislocation glide can be simulated by direct molecular dynamics (MD) simulations³⁻⁵ (as illustrated in Fig. 1 (a) and left track in Fig. 1 (b)),

but only at high temperatures and/or high strain rates. To access more general conditions, one can use the Transition State Theory (TST)⁶, which expresses the average forward dislocation velocity $v$ at a temperature $T$ under an applied stress tensor $[σ]$:

$$v = vL \exp \left( - \frac{\Delta G([σ], T)}{k_B T} \right), \qquad (1)$$

with $v$ an attempt frequency, $L$ the dislocation length, $k_B T$ the thermal energy and $ΔG$ the Gibbs activation energy for kink-pair nucleation. In BCC metals, the activation free enthalpy, $ΔG$, depends a priori on all components and sign of the applied stress tensor due to non-Schmid

<footnote n="1">¹ASNR/PSN-RES/SEMIA/LSMA Centre d'études de Cadarache, F-13115 Saint Paul-lez-Durance, France. ²Univ. Lyon, UCBL, Institut Lumière Matière, UMR CNRS 5306, F-69622 Villeurbanne, France. ³Univ. Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, SRMP, F-91191 Gif-sur-Yvette, France. ⁴Univ. Lyon, INSA Lyon, UCBL, MATEIS, UMR CNRS 5510, F-69621 Villeurbanne, France. ⁵Aix-Marseille Univ., CINaM, UMR CNRS 7325, Campus de Luminy, F-13288 Marseille, France. ⁶Department of Mechanical Engineering, University of Michigan, Ann Arbor 48109 Michigan, USA. ⁷Max Planck Institute for Sustainable Materials, Max-Planck-Straße 1, 40237 Düsseldorf, Germany. ✉ e-mail: arnaud.allera@asnr.fr; david.rodney@univ-lyon1.fr</footnote>
effects<sup>7,8</sup>. However, in this work, we will only focus on the effect of a resolved shear stress applied in the dislocation glide plane.

Due to the difficulty in computing Gibbs activation energies, ad hoc simplifications are often used. The most classical is the harmonic TST (HTST)<sup>10</sup> (see central track in Fig. 1 (b)). Both the activation entropy $\Delta S_h$ and enthalpy $\Delta H_h$ are then independent of $T$. $\Delta S_h$ is computed by diagonalizing the Hessian matrix of the system in its zero-Kelvin initial state (a straight screw dislocation in a Peierls valley) and activated state (an unstable double-kinked dislocation in-between Peierls valleys, illustrated in Fig. 1a)<sup>11</sup>. However, the diagonalizations are very demanding in computational resources for systems containing the $N \sim 10^2$ atoms needed to model dislocations<sup>11</sup>. As a result, further approximations are commonly applied. The most common is to simply neglect the stress dependence of $\Delta S_h$, which becomes a constant that can be absorbed into the attempt frequency $\nu' = \nu \exp(\Delta S_h/k_B)$. However, there is evidence that activation entropies in dislocation-related processes may exhibit a stress dependence, including dislocation nucleation<sup>12,13</sup>, obstacle bypass<sup>14</sup>, cross-slip<sup>15</sup>, and the glide of straight dislocations<sup>16</sup>. It is important to note here that these calculations were performed using empirical interatomic potentials. An alternative approach to incorporate stress dependence without explicitly computing $\Delta S$ is to invoke the Meyer-Neldel (MN) law, or compensation rule<sup>17</sup>, which states that the activation entropy of a process is proportional to its activation enthalpy, $\Delta S = \Delta H/T_{MN}$, where $T_{MN}$ is a characteristic temperature. Despite its widespread use, the MN law lacks a firm physical foundation, and counterexamples have been reported in both experiments<sup>18</sup> and computational simulations<sup>19–22</sup>. Here also, it should be emphasized that these computational studies were based on empirical potentials.

Empirical potentials, such as embedded atom method (EAM) potentials, have significantly advanced the modeling of BCC screw dislocations. However, most fail to accurately reproduce the potential energy landscape of screw dislocations obtained from density functional theory (DFT) calculations<sup>7</sup> and the temperature dependence of elastic constants compared to experimental data<sup>23</sup>. Machine-learning interatomic potentials (MLIPs) provide a promising alternative to overcome these limitations<sup>24</sup> at a computational cost that has now become manageable. However, while some recently developed MLIPs for Fe have been explicitly trained on dislocation configurations<sup>25–27</sup>, their datasets remain limited and often focus primarily on zero-Kelvin properties.

In this work, we investigate the influence of vibrational entropy on dislocation mobility without relying on a priori simplifications, using state-of-the-art interatomic potentials. We modeled isolated dislocations in BCC metals and directly computed their (possibly anharmonic) Gibbs activation energy for kink-pair nucleation using the linear-scaling projected average force integrator (PAFI)<sup>28</sup>. This offers an alternative computational route to estimate the flow stress, following the path highlighted in Fig. 1b. To accurately capture activation thermodynamics, we developed MLIPs for two common BCC metals, Fe and W, by refining the models from ref. 26 and expanding the training dataset to include a wider variety of dislocation configurations. Further details on the simulation setup, computational methods, and MLIPs are provided in the Methods section.

# Results

## Machine-learning interatomic potentials

We begin by examining the Gibbs activation energy for kink-pair nucleation in Fe. To this end, we developed a MLIP for Fe, building upon the framework established in ref. 26 and successfully applied to free energy calculations of defects in BCC metals<sup>23,29</sup>. The training dataset was significantly expanded to include a wide range of dislocation configurations. Notably, in addition to the Peierls barrier, the dataset incorporates the hard-to-split path<sup>7</sup> and single-kink configurations<sup>30</sup>, which are considered for the first time. These

a

<figure><img src="image_4.png" /></figure>

b

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Overview of dislocation glide simulations and associated computational methods.** a A screw dislocation migrating from one Peierls valley to the next by a kink-pair mechanism. The 3D atomic structure was rendered with Ovito<sup>9</sup>. The upper part of the crystal is not shown for clarity. b Overview of different computational approaches to predict the flow stress.</figcaption></figure>

configurations were evaluated both at 0 K and finite temperatures, with non-equilibrium states sampled using constrained MD simulations (see Supplementary Material for details). The resulting MLIP accurately reproduces the Peierls mechanism, core trajectory, and kink-pair nucleation enthalpy (see Fig. 2 (a-c)), as well as core eigen-strains (see Supplementary Material), demonstrating excellent agreement with DFT data. For comparison, we included in Fig. 2 predictions of a classical EAM potential for Fe<sup>11</sup>, which will be further used in “EAM potentials”.

Using the MLIP, we employed the PAFI method to compute minimum Gibbs energy paths (MGEP) without any a priori assumptions. The resulting Gibbs activation energies, obtained at 300 and 600 MPa, are shown in Fig. 2d. At low temperatures, the enthalpy decreases linearly with increasing temperature, indicating a harmonic regime that extends up to $T_0 \approx 100$ K. Beyond this threshold, the Gibbs activation energy deviates from the harmonic extrapolation and decreases more gradually. The discrepancy between HTST and anharmonic calculations is further illustrated in Fig. 2e and discussed in Supplementary Material. This difference primarily arises from a change in the Gibbs energy profile: above $T_0$, the position of the barrier minimum shifts from $z=0$ to a finite coordinate $z_1$, while the barrier maximum remains mostly fixed at a coordinate $z_2$ (see Fig. 2e). Unfortunately, we could not identify any apparent structural changes that may account for the lowering of the Gibbs free energy in the finite-temperature configurations generated by PAFI. In Supplementary Material, we detail a statistical outlier analysis of the shifted transition pathway sampled by PAFI, demonstrating that the sampled configurations remain in the interpolation domain of the MLIP.

To test the validity of the harmonic assumption along the transition path and extend the harmonic regime, we use the variational HTST
<figure><img src="image_6.png" /><figcaption>**Fig. 2 |** MLIP predicts a harmonic behavior and a stress-independent activation entropy. **a** Peierls barrier for a 1/2[111] straight screw dislocation and **(b)** corresponding core trajectory, computed with the present MLIP for Fe and a classical EAM potential¹¹, compared to DFT data³¹. **c** Kink-pair nucleation enthalpy computed with the EAM and MLIP, compared to a line tension (LT) model parameterized on DFT calculations³². **d** Gibbs activation energy for kink-pair nucleation computed with the MLIP. PAFI calculations are compared to variational and non-</figcaption></figure>

variational HTST predictions. Shaded regions represent the error estimation of the PAFI method, as described in Supplementary Material. **e** Examples of Gibbs energy profiles at 600 MPa showing a shift of the minimum Gibbs energy configuration. **f** Dislocation velocity computed by direct MD simulations with the MLIP (circles) under a constant applied shear stress, compared to Eq. (1) parameterized with data from (c-d) with $v = 3.8 \times 10^9$ Hz for HTST (dotted lines) and $9.2 \times 10^{10}$ Hz for variational HTST (dashed lines).

(VHTST)³³, where the Gibbs energy barrier is computed harmonically between the configurations that maximize the barrier height, a method recently employed to model dislocation nucleation kinetics³⁴. We thus approximate $\Delta G_{VHTST} = \Delta H_h(z_2) - T\Delta S_h(z_2)$ if $T < T_o$ and $\Delta G_{VHTST} = \Delta H_h(z_2) - \Delta H_h(z_1) - T(\Delta S_h(z_2) - \Delta S_h(z_1))$ if $T > T_o$, where $\Delta H_h(z)$ and $\Delta S_h(z)$ are the OK enthalpy and entropy differences between the configuration at reaction coordinate $z$ and the initial configuration (at $z=0$). The result is reported in Fig. 2d as dotted lines, showing a very good agreement with the anharmonic PAFI calculations over the entire temperature range considered here.

Perhaps the most striking result is the remarkably weak dependence of the activation entropy on the applied stress: at low temperature, in the harmonic regime, $\Delta S(z_2) = 6.3 \text{ kB}$, while above $T_o$, $\Delta S(z_2) - \Delta S(z_1) = 1.6 \text{ kB}$. In the harmonic regime, the activation entropy is thus independent of both temperature and applied stress, providing, for the first time, a direct justification for the commonly used simplification of a constant entropy.

To further validate our calculations, we conducted MD simulations with the MLIP for different resolved shear stresses and temperatures. The measured average dislocation velocities are presented in Fig. 2 (f) and compared to the predictions of both HTST, $v = \nu L \exp(-\Delta G_h/k_B T)$, and VHTST, $v = \nu L \exp(-\Delta G_{VHTST}/k_B T)$, where $v$ is the only fitting parameter. Both approximations yield nearly identical velocity predictions, as differences in entropy are absorbed into the fitted prefactor $v$, while the enthalpic change is negligible, given that $\Delta H_h(z_1) \ll \Delta H_h(z_2)$. The strong agreement between MD and theoretical predictions confirms the harmonic nature of the glide transition and the absence of significant stress dependence of the activation entropy for this MLIP.

To evaluate the generality of the activation enthalpy and entropy behavior observed in Fe, we developed a second MLIP for W. As in Fe, a database built in previous works²³,²⁶,³⁵ was enriched with finite-temperature dislocation configurations. The corresponding results, presented in Supplementary Material, confirm the existence of a harmonic regime up to at least 100 K with an activation entropy independent of the applied stress. Remarkably, despite the significant differences between Fe and W -such as their elastic constants - the harmonic activation entropy in W is found to be approximately $8k_B$, closely matching the value obtained for Fe. This unexpected similarity

suggests a degree of universality in the contribution of vibrational entropy to kink pair nucleation across BCC metals.

## EAM potentials

The strikingly constant entropy obtained with the MLIPs for Fe and W contrasts with the results obtained in the literature using empirical potentials and mentioned in the Introduction. To verify this point, we used two EAM potentials, one for Fe¹¹, the other for W³⁶, widely used to simulate dislocation glide at the atomic scale. We computed their Gibbs activation energies for kink-pair nucleation, and since the calculations were much faster than with the MLIPs, we could consider a wide range of applied stresses. Figure 3 considers the case of Fe, while the results in W are shown in Supplementary Material. Figure 3a shows the Gibbs activation energies at 200 and 500 MPa. We see that $\Delta G$ exhibits a temperature dependence more complex than with the MLIP, with a drastic departure from the harmonic prediction (dashed lines) above temperatures as low as 20 K. Note that we have checked that for the moderate stresses and temperatures considered here, classically discussed anharmonic effects, temperature dependence of the elastic moduli¹³ and thermal dilatation²² do not affect the energy barriers (see Supplementary Material). The main effect is rather a widening of the activated kink pair²⁸, which reflects a deviation of the MGEPs at finite temperature, as detailed in Supplementary Material.

Although the HTST approximation appears inapplicable beyond about 20 K with the present EAM potential, the evolution of $\Delta G$ shown in Fig. 3b can still be approximated by a linear relation with reasonable accuracy below 100 K (dashed lines in Fig. 3b). This defines an effective temperature-independent entropy, $\Delta S_{eff}(\tau)$. Interestingly, we see in Fig. 3 (b) that the initial slope of $\Delta G$ becomes steeper at higher stresses, which means that $\Delta S_{eff}(\tau)$ increases with the applied stress. Plotted against $\Delta H_h(\tau)$ in Fig. 3c, we find a linear variation, which follows an inverse Meyer-Neldel law, or reinforcement effect: $\Delta S_{eff}$ decreases with increasing $\Delta H_h$, instead of increasing as assumed by the MN compensation law. Note also that $\Delta S_{eff}$ varies by about $10k_B$ between 0 and 700 MPa, in stark contrast with the results obtained with the MLIPs. The emergence of an inverse Meyer-Neldel law is consistent with recent works presented in refs. 21,22, indicating that it might be a general trend of EAM potentials. This conclusion is reinforced by calculations
<figure><img src="image_7.png" /><figcaption>**Fig. 3 | Empirical potential predicts strong anharmonicity and an inverse Meyer-Neldel relation.** a Gibbs activation energy as a function of temperature at 200 and 500 MPa, obtained with PAFI (circles and full lines) and the HTST approximation using Hessian matrix diagonalizations (dashed lines). Shaded regions represent the error estimation of the PAFI method, as described in Supplementary Material. Lines are colored based on the associated stress value.</figcaption></figure>

b Gibbs activation energy as a function of temperature for a large range of stresses and temperatures. Linear fits up to 100 K with slope –$\Delta S_{\text{eff}}$ are shown as dashed lines. Note that calculations presented in (a) were performed with a high number of independent samplings (ranging from 150–250), in order to accurately capture subtle variations of Gibbs activation energy at very low temperature. c Effective entropy $\Delta S_{\text{eff}}$ as a function of enthalpy.

performed with the EAM potential for W. As presented in Supplementary Material, the potential also predicts a harmonic regime limited to very low temperatures with again an inverse MN rule between the enthalpy and a strongly varying entropy.

## From Gibbs energy to yield stress

To illustrate the direct impact of entropic effects on macroscopic properties, we employed the analytical model from ref. 37 to predict the thermally activated yield stress in Fe and W. The predictions are compared with experimental data from refs. 38,39 for Fe and ref. 40 for W. The model requires no fitting parameters apart from the attempt frequency and can be fully parameterized using either MLIPs or EAM potentials (see the Methods section for details). For calculations using EAM potentials, we incorporated the inverse MN relation, $\Delta S = - \Delta H/T_{\text{MN}}$, with characteristic temperatures of 406 K for Fe and 1078 K for W as obtained in the PAFI calculations “EAM potentials”. In contrast, for calculations using MLIPs, and in agreement with the results presented in “Machine-learning interatomic potentials”, we assumed a constant activation entropy of $\Delta S = 6.3 k_B$ for both metals. To apply the model under conditions representative of experimental settings, one parameter difficult to evaluate is the density of mobile dislocations. To address this, we considered a range of densities from $10^6$ to $10^{11} \text{ m}^{-2}$, using $10^8 \text{ m}^{-2}$ as a reference value, in agreement with estimates from the literature$^{41–43}$.

The results are presented in Fig. 4. As with all previous EAM and DFT calculations, the present MLIPs overestimate the Peierls stress at 0 K compared to experiments. To eliminate this well-known discrepancy$^{11,44,45}$, all curves were normalized by the extrapolated zero-Kelvin yield stress. Both MLIP predictions exhibit excellent agreement with experimental data. Notably, they accurately capture the athermal temperature, above which plasticity ceases to be thermally activated. For the EAM potentials, it is interesting to observe that the prediction for Fe is rather accurate, although with deviations at both low and high temperatures, whereas in W, the EAM potential significantly underestimates the temperature dependence of the yield stress. We believe that the agreement of the EAM potential in Fe results from a compensation effect: the kink-pair nucleation enthalpy is underestimated, as seen in Fig. 2, while the entropy is overestimated, as shown in Fig. 3. In W, the enthalpy is better reproduced, but the entropy remains overestimated, leading to a slower decrease in the yield stress. Regardless, these results clearly highlight the ability of MLIPs to predict macroscopic yield stress without the need for fitting parameters.

## Discussion

The main conclusion of this work may seem unremarkable at first: MLIPs reveal that the activation entropy for dislocation glide in BCC Fe and W is constant, as assumed in the simplest models. However, this finding contrasts sharply with previous studies based on EAM potentials$^{12–16,21,34}$, suggesting a strongly varying activation entropy indicative of pronounced anharmonic effects. The present calculations indicate that this variability is likely an artifact of EAM potentials, which tend to produce a more rugged energy landscape compared to MLIPs and DFT. Nevertheless, we emphasize that potential energy landscape artifacts may still arise even with MLIPs$^{46}$, and care must be taken to validate their reliability in each application.

Vibrational entropy arises from subtle variations in the vibrational density of states, which EAM potentials struggle to capture due to their simplified formulation. In contrast, MLIPs exhibit a consistently harmonic behavior, as demonstrated by the accurate entropy predictions across the entire temperature range when using a variational HTST approach. This stands in stark contrast to the pronounced anharmonicity inherent in EAM potentials.

While our results highlight several limitations of classical EAM potentials compared to MLIPs, we acknowledge that a formal demonstration of the origin of these limitations remains elusive. What is clear, however, is that the central-force nature of EAM potentials leads to enhanced shear softening. Additional artifacts have also been reported, such as a non-monotonic temperature dependence of elastic constants$^{23}$ and an early departure from harmonicity$^{28}$. Although some rugged features in energy landscapes, such as those reported in recent studies on metallic glasses$^{47,48}$, may correspond to real physical phenomena, our findings highlight the need for caution in interpreting such features, as they may also stem from artifacts due to the potential. This underscores the importance of validating potential energy landscapes against DFT data or using well-trained MLIPs, and calls for more systematic cross-comparisons between different potential formalisms in future works.

As a final note, we recall that early theories on entropic effects$^{49–51}$ proposed a connection between entropy and the temperature dependence of macroscopic properties, particularly the shear modulus $\mu$. Zener$^{49}$, for instance, suggested that $\Delta S \propto d\mu/dT$. To investigate this hypothesis, we computed the temperature dependence of the shear modulus using both the EAM potential and the MLIP for Fe (see Supplementary Material). While the EAM potential predicts an anomalous increase in shear modulus with temperature-potentially aligning with an inverse compensation effect in Zener’s theory-the
<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Yield stress model parameterized on atomistic calculations, compared to experimental data.** The model from ref. 37 was parameterized using EAM potentials or MLIPs, in Fe and W. Experimental data are reproduced from refs. 38,39</figcaption></figure>

for Fe and ref. 40 in W. Shaded areas indicate the variation range due to changes in dislocation density.

MLIP predicts a decreasing shear modulus, a trend inconsistent with the absence of entropic effects in this potential. These findings thus suggest that inferring an enthalpy-entropy relationship from macroscopic property variations is challenging. Instead, entropy changes likely emerge from complex changes in vibrational spectra, as previously observed in amorphous solids²⁰.

Although this study focuses on metal plasticity, the results are expected to extend to other crystalline defects, such as point defects and grain boundaries, whose kinetics may also be influenced by entropic effects. Caution must therefore be exercised when analyzing finite-temperature kinetics predicted from the atomic scale.

## Methods

### Machine-learning interatomic potentials

The machine-learning interatomic potentials are based on the Quadratic Noise Machine Learning (QNML) formalism introduced in refs. 23,26,29. This approach provides a balanced compromise between computational efficiency and predictive accuracy. Local atomic environments are represented by bispectrum SO(4) descriptors with maximum angular momentum $J_{max} = 4$ and dimension $K = 55$. The cutoff distance is set to 4.7 Å for Fe and 5.3 Å for W. The fitting of the potentials was carried out using the MILADY package⁵². We used weights in the objective loss function to control accuracy in different subsets of the database and optimize properties of interest. We extended existing general-purpose fitting databases for Fe and W²³,²⁶ with dislocation configurations. For consistency, we ensured that the new configurations were computed with the same DFT parameterization as the original databases, using the VASP code. The Mahalanobis

outlier distance is used in Supplementary Material to detect out-of-distribution configurations⁵³.

### Atomistic simulations of screw dislocations

The simulation cell contains a screw dislocation of length $L = 40$ b, where b is the Burgers vector of the dislocation, defined as $b = \frac{1}{2}[111]$. The dislocation glides in a {110} central plane with a glide distance $D = 50$ b and a cell height perpendicular to the glide plane of $H = 30$ b. Periodic boundary conditions were applied in the glide plane to form a periodic array of dislocations⁵⁴. The total number of atoms in the simulation cell is 96,000.

Stress-controlled molecular dynamics simulations were performed using LAMMPS⁵⁵ with the MLIPs. The simulations lasted from 0.2 to 0.6 ns to ensure the computation of well-converged average dislocation velocities. The simulation time step was set to 1 fs, with a total computational cost between $3 \times 10^4$ and $9 \times 10^4$ CPU hours per condition.

### Computation of Gibbs energy barriers

Gibbs energy barriers were calculated using the linear-scaling Projected Average Force Integrator (PAFI) method²⁸. Transition pathways were first relaxed using the nudged elastic band (NEB) method⁵⁶ to a force tolerance of $10^{-3}$ eV Å⁻¹. A custom parallel nudging force, designed to promote equidistant energy differences along the path, was implemented in LAMMPS⁵⁵.

The PAFI approach enables the calculation of the anharmonic Gibbs activation energy at $O(N)$ cost for kink-pair nucleation in bcc metals under finite applied stress along the 0K reaction coordinate,
without requiring any a priori assumptions about finite temperature
changes of the pathways. Anharmonic Gibbs energy calculations were
performed using the PAFI package<sup>28,57,58</sup>. Within this method, the
configurations at a given applied stress and specific reaction coordi-
nate are sampled at all temperatures from the same hyperplane, which
is perpendicular to the initial NEB path. Furthermore, after sampling,
configurations are relaxed in their hyperplanes to eliminate trajec-
tories that might have escaped their initial potential energy valley. The
computational cost ranged from 5 × 10<sup>4</sup> to 2.5 × 10<sup>5</sup> CPU hours per
condition, mainly depending on the number of independent sam-
plings, which typically ranges from 5 to 40 and was optimized to
control numerical errors (see Supplementary Material for a discussion
of error estimation in PAFI calculations).

Hessian matrix diagonalizations were used to compute the har-
monic activation entropy within the framework of the harmonic
transition state theory. The Hessian matrices were generated using
LAMMPS and subsequently diagonalized with the Phony code to
evaluate the phonon spectrum<sup>59–61</sup>. The total computational cost for
these calculations amounted to approximately 5 × 10<sup>4</sup> CPU-hours per
atomic system using the MLIP.

Analytical yield stress model

The thermally activated yield stress in Fe and W was predicted using
the analytical model developed in ref. 37, which expresses the nor-
malized yield stress as:

$$
\frac{\sigma_Y(T)}{\sigma_Y(T=0K)} = \left( 1 - \left( \frac{k_B T}{\Delta H_0 \left(1 - \frac{T}{T_{MN}}\right)} \ln \left( \frac{\dot{\epsilon}_p}{\sqrt{\rho_D \nu_D \lambda_p}} \right) \right)^{\frac{1}{q}} \right)^{\frac{1}{p}} . \quad (2)
$$

This expression derives from a Kocks representation of the acti-
vation enthalpy:

$$
\Delta H = \Delta H_0 \left( 1 - \left( \frac{\sigma_Y(T)}{\sigma_Y(T=0K)} \right)^p \right)^q \quad (3)
$$

and a MN law for the activation entropy:

$$
\Delta S = \frac{\Delta H}{T_{MN}}. \qquad (4)
$$

Activation enthalpies were fitted to NEB calculations using the different potentials. For Fe, we used ΔH₀ = 0.6 eV, p = 0.53 and q = 1.06 for the EAM potential, and ΔH₀ = 0.8 eV, p = 0.87 and q = 1.33 for the MLIP. For W, we used ΔH₀ = 1.6 eV, p = 0.68 and q = 1.02 for the EAM potential, and ΔH₀ = 1.54 eV, p = 0.86 and q = 1.43 for the MLIP. The harmonic behavior predicted by the MLIPs was obtained in the limit T_MN → ∞ while the inverse MN law predicted by the EAM potentials was represented with T_MN = -406 K in Fe and -1078 K in W. The strain rates were set to the experimental values, ˙εₚ = 5.6 × 10⁻⁴ s⁻¹ in Fe, and 8.5 × 10⁻⁴ s⁻¹ in W. The distance between Peierls valleys is λₔ = 2.31 Å in Fe, 2.6 Å in W. The attempt frequency was set to νₔ = 6.15 × 10¹³ s⁻¹ in Fe and 1 × 10¹³ s⁻¹ in W.

Reporting summary

Further information on research design is available in the Nature
Portfolio Reporting Summary linked to this article.

Data availability

The data generated in this study have been deposited in the Zenodo
repository 10.5281/zenodo.15746898<sup>62</sup>. The repository includes the Fe
and W MLIPs and a minimal dataset to reproduce the results. Source
data are provided with this paper.

Code availability

The MILADY package is open-source software distributed under the
Apache Software License (ASL)⁵². The PAFI package is available on
GitHub⁵⁷ and is also included in the EXTRA-FIX package of the LAMMPS
software.

**References**

1. Caillard, D. & Martin, J. L. *Thermally Activated Mechanisms In Crystal Plasticity* 1st edn, Vol. 8 (Elsevier, 2003).

2. Po, G. et al. A phenomenological dislocation mobility law for bcc metals. *Acta Mater.* **119**, 123–135 (2016).

3. Domain, C. & Monnet, G. Simulation of screw dislocation motion in iron by molecular dynamics simulations. *Phys. Rev. Lett.* **95**, 215506 (2005).

4. Chaussidon, J., Fivel, M. & Rodney, D. The glide of screw disloca-
tions in bcc Fe: atomistic static and dynamic simulations. *Acta
Mater.* **54**, 3407–3416 (2006).

5. Zepeda-Ruiz, L. A., Stukowski, A., Oppelstrup, T. & Bulatov, V. V. Probing the limits of metal plasticity with molecular dynamics simulations. *Nature* **550**, 492–495 (2017).

6. Laidler, K. J. & King, M. C. Development of transition-state theory. *J. Phys. Chem.* **87**, 2657–2664 (1983).

7. Dezerald, L., Rodney, D., Clouet, E., Ventelon, L. & Willaime, F. Plastic anisotropy and dislocation trajectory in BCC metals. *Nat. Commun.* **7**, 1–7 (2016).

8. Kraych, A. et al. Non-glide effects and dislocation core fields in BCC metals. *NJP Comput. Mater.* **5**, 1–8 (2019).

9. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO—the open visualization tool. *Model. Simul. Mater. Sci. Eng.* **18**, 015012 (2009).

10. Vineyard, G. H. Frequency factors and isotope effects in solid state rate processes. *J. Phys. Chem. Solids* **3**, 121–127 (1957).

11. Proville, L., Rodney, D. & Marinica, M.-C. Quantum effect on ther-
mally activated glide of dislocations. *Nat. Mater.* **11**, 845–849 (2012).

12. Warner, D. H. & Curtin, W. Origins and implications of temperature-
dependent activation energy barriers for dislocation nucleation in
face-centered cubic metals. *Acta Mater.* **57**, 4267–4277 (2009).

13. Ryu, S., Kang, K. & Cai, W. Entropic effect on the rate of dislocation nucleation. *Proc. Natl. Acad. Sci. USA* **108**, 5174–5178 (2011).

14. Saroukhani, S., Nguyen, L., Leung, K., Singh, C. & Warner, D. Har-
nessing atomistic simulations to predict the rate at which disloca-
tions overcome obstacles. *J. Mech. Phys. Solids* **90**, 203–214 (2016).

15. Esteban-Manzanares, G., Santos-Güemes, R., Papadimitriou, I., Martínez, E. & LLorca, J. Influence of the stress state on the cross-slip free energy barrier in Al: an atomistic investigation. *Acta Mater.* **184**, 109–119 (2020).

16. Gilbert, M., Schuck, P., Sadigh, B. & Marian, J. Free energy gen-
eralization of the Peierls potential in iron. *Phys. Rev. Lett.* **111**,
095502 (2013).

17. Yelon, A., Movaghar, B. & Crandall, R. Multi-excitation entropy: its role in thermodynamics and kinetics. *Rep. Prog. Phys.* **69**, 1145 (2006).

18. Meiling, H. & Schropp, R. The inverse Meyer–Neldel rule in thin-film transistors with intrinsic heterogeneous silicon. *Appl. Phys. Lett.* **74**, 1012–1014 (1999).

19. Koziatek, P., Barrat, J.-L., Derlet, P. & Rodney, D. Inverse Meyer-Neldel behavior for activated processes in model glasses. *Phys. Rev. B* **87**, 224105 (2013).

20. Gelin, S., Champagne-Ruel, A. & Mousseau, N. Enthalpy-entropy compensation of atomic diffusion originates from softening of low frequency phonons. *Nat. Commun.* **11**, 3977 (2020).

21. Zotov, N. & Grabowski, B. Entropy of kink pair formation on screw dislocations: an accelerated molecular dynamics study. *Model. Simul. Mater. Sci. Eng.* **30**, 065004 (2022).
22. Wang, Y. & Cai, W. Stress-dependent activation entropy in thermally activated cross-slip of dislocations. *Proc. Natl. Acad. Sci. USA* **120**, e2222039120 (2023).

23. Zhong, A. et al. Anharmonic thermo-elasticity of tungsten from accelerated Bayesian adaptive biasing force calculations with data-driven force fields. *Phys. Rev. Mater.* **7**, 023802 (2023).

24. Deringer, V. L., Caro, M. A. & Csányi, G. Machine learning interatomic potentials as emerging tools for materials science. *Adv. Mater.* **31**, 1902765 (2019).

25. Byggmästar, J. et al. Multiscale machine-learning interatomic potentials for ferromagnetic and liquid iron. *J. Phys. Condens. Matter* **34**, 305402 (2022).

26. Goryaeva, A. M. et al. Efficient and transferable machine learning potentials for the simulation of crystal defects in bcc Fe and W. *Phys. Rev. Mater.* **5**, 103803 (2021).

27. Meng, F.-S. et al. General-purpose neural network interatomic potential for the α-iron and hydrogen binary system: toward atomic-scale understanding of hydrogen embrittlement. *Phys. Rev. Mater.* **5**, 113606 (2021).

28. Swinburne, T. D. & Marinica, M.-C. Unsupervised calculation of free energy barriers in large crystalline systems. *Phys. Rev. Lett.* **120**, 135503 (2018).

29. Zhong, A. et al. Unraveling temperature-induced vacancy clustering in tungsten: from direct microscopy to atomistic insights via data-driven Bayesian sampling. *PRX Energy* **4**, 013008 (2025).

30. Ventelon, L., Willaime, F. & Leyronnas, P. Atomistic simulation of single kinks of screw dislocations in α-fe. *J. Nucl. Mater.* **386**, 26–29 (2009).

31. Bienvenu, B., Dezerald, L., Rodney, D. & Clouet, E. Ab initio informed yield criterion across body-centered cubic transition metals. *Acta Mater.* **236**, 118098 (2022).

32. Dezerald, L., Proville, L., Ventelon, L., Willaime, F. & Rodney, D. 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 (2015).

33. Truhlar, D. G. & Garrett, B. C. Variational transition state theory. *Annu. Rev. Phys. Chem.* **35**, 159–189 (1984).

34. Bagchi, S. & Perez, D. Anomalous entropy-driven kinetics of dislocation nucleation. *Nat. Commun.* **16**, 912 (2025).

35. Zhong, A. *Machine Learning And Adaptive Sampling To Predict Finite-Temperature Properties In Metallic Materials At The Atomic Scale*. https://theses.hal.science/tel-04764694v1 (2024).

36. Marinica, M.-C. et al. Interatomic potentials for modelling radiation defects and dislocations in tungsten. *J. Phys. Condens. Matter* **25**, 395502 (2013).

37. Clouet, E., Bienvenu, B., Dezerald, L. & Rodney, D. Screw dislocations in bcc transition metals: from ab initio modeling to yield criterion. *Comptes Rendus. Phys.* **22**, 83–116 (2021).

38. Keh, A. & Nakada, Y. Plasticity of iron single crystals. *Can. J. Phys.* **45**, 1101–1120 (1967).

39. Zwiesele, S. & Diehl, J. Temperature and strain rate dependence of the macro yield stress of high purity iron single crystals, in *Strength of Metals and Alloys*, (eds. Haasen, P., Gerold, V. & Kostorz, G.) 59–64 (Elsevier, 1979).

40. Brunner, D. & Glebovsky, V. The plastic properties of high-purity w single crystals. *Mater. Lett.* **42**, 290–296 (2000).

41. Brunner, D. & Glebovsky, V. Analysis of flow-stress measurements of high-purity tungsten single crystals. *Mater. Lett.* **44**, 144–152 (2000).

42. Brunner, D. Comparison of flow-stress measurements on high-purity tungsten single crystals with the kink-pair theory. *Mater. Trans.*, **JIM 41**, 152–160 (2000).

43. Caillard, D. Kinetics of dislocations in pure Fe. part ii. in situ straining experiments at low temperature. *Acta Mater.* **58**, 3504–3515 (2010).

44. Freitas, R., Asta, M. & Bulatov, V. V. Quantum effects on dislocation motion from ring-polymer molecular dynamics. *NPJ Comput. Mater.* **4**, 55 (2018).

45. Proville, L., Ventelon, L. & Rodney, D. Prediction of the kink-pair formation enthalpy on screw dislocations in α-iron by a line tension model parametrized on empirical potentials and first-principles calculations. *Phys. Rev. B* **87**, 144106 (2013).

46. Deng, B. et al. Systematic softening in universal machine learning interatomic potentials. *NPJ Comput. Mater.* **11**, 1–9 (2025).

47. Zella, L., Moon, J. & Egami, T. Ripples in the bottom of the potential energy landscape of metallic glass. *Nat. Commun.* **15**, 1358 (2024).

48. Li, H., Xiao, H., Egami, T. & Fan, Y. Infinitely rugged intra-cage potential energy landscape in metallic glasses caused by many-body interaction. *Mater. Today Phys.* **49**, 101582 (2024).

49. Zener, C. Theory of $D_0$ for atomic diffusion in metals. *J. Appl. Phys.* **22**, 372–375 (1951).

50. Schoeck, G. The activation energy of dislocation movement. *Phys. Status Solidi (b)* **8**, 499–507 (1965).

51. Kocks, U. F., Argon, A. S. & Ashby, M. F. Thermodynamics and kinetics of slip, In *Progress in Materials Science*, (eds. Chalmers, J. W., Christian, T. B., Massalski) 1–281 (Oxford, 1975).

52. Milady. What is MILADY?. https://ai-atoms.github.io/milady/. (2022).

53. Goryaeva, A. M. et al. Reinforcing materials modelling by encoding the structures of defects in crystalline solids into distortion scores. *Nat. Commun.* **11**, 4691 (2020).

54. Bacon, D., Osetsky, Y. & Rodney, D. Chapter 88 dislocation-obstacle interactions at the atomic level. In *Chapter 88 Dislocation-Obstacle Interactions at the Atomic Level, Vol. 15 of Dislocations in Solids* (eds. Hirth, J. & Kubin L.) 1–90 (Elsevier, 2009).

55. Thompson, A. P. et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comp. Phys. Comm.* **271**, 108171 (2022).

56. Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. *J. Chem. Phys.* **113**, 9901–9904 (2000).

57. Swinburne, T. D., Allera, A., McMarinica. tomswinburne/pafi: 2025 v0.9.9. Zenodo https://doi.org/10.5281/zenodo.15730212 (2025).

58. Sato, Y., Swinburne, T., Ogata, S. & Rodney, D. Anharmonic effect on the thermally activated migration of {1012} twin interfaces in magnesium. *Mater. Res. Lett.* **9**, 231–238 (2021).

59. Marinica, M. C. & Willaime, F. Orientation of interstitials in clusters in α-fe: a comparison between empirical potentials. *Solid State Phenom.* **129**, 67–74 (2007).

60. Soulié, A., Bruneval, F., Marinica, M.-C., Murphy, S. & Crocombette, J.-P. Influence of vibrational entropy on the concentrations of oxygen interstitial clusters and uranium vacancies in nonstoichiometric $UO_2$. *Phys. Rev. Mater.* **2**, 083607 (2018).

61. Berthier, F. et al. Order-disorder or phase-separation transition: analysis of the Au-Pd system by the effective site energy model. *Phys. Rev. B* **99**, 014108 (2019).

62. Allera, A. et al. Supporting files for “activation entropy of dislocation glide in body-centered cubic metals from atomistic simulations”. Zenodo. https://doi.org/10.5281/zenodo.15746898 (2025).

## Acknowledgements

Emmanuel Clouet, Lisa Ventelon and Thomas Leveau are gratefully acknowledged for insightful discussions. AA acknowledges financial support from the Cross-Disciplinary Program on Numerical Simulation of the French Alternative Energies and Atomic Energy Commission (CEA). A.A., A.M.G., T.D.S and M.-C.M. acknowledge support from GENCI—(Jean-Zay/CINES/CCRT) computer centre under Grant No. A0170906973. T.D.S gratefully acknowledges support from ANR grant ANR-23-CE46-0006-1, IDRIS allocations A0120913455, and Euratom Grant No. 633053.
**Author contributions**

A.A., T.D.S., M.-C.M. and D.R. conceived the study. B.B., M.-C.M., A.M.G. and A.A. created the training dataset and performed DFT calculations. A.A., M.-C.M. and A.M.G. developed and fitted the machine learning potentials. A.A. and T.D.S. conducted the MD and PAFI simulations. A.A. and M.-C.M. performed harmonic entropy calculations. A.A and B.B. adapted the yield model. All authors contributed to the interpretation of the results and participated in scientific discussions throughout the study. A.A. and D.R. drafted the initial manuscript and subsequent revisions, which were reviewed by all authors.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-62390-w.

**Correspondence** and requests for materials should be addressed to Arnaud Allera or David Rodney.

**Peer review information** *Nature Communications* thanks Yifan Wang, Xinran Zhou and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. [A peer review file is available].

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025