As a baseline for subsequent comparison, first we compute the Peierls stress using fully classical simulations. Defined as the minimal stress required to make a dislocation move in the limit of zero temperature, Peierls stress can be computed in a static simulation in which stress or strain is increased in increments, each increment followed by a full relaxation of atom positions. The Peierls stress is taken to be just one increment below the stress at which the dislocation begins to move. Our so-computed “static” Peierls stress is τ P = (0.98 ± 0.01) GPa for a \({\textstyle{1 \over 2}}\left\langle {111} \right\rangle\) screw dislocation moving on a {112} plane in the so-called twinning direction.24 This value is close to the ones previously reported for the same potential10,11 (≈0.91 GPa). For reasons explained in Supplementary Information 1, the RPMD method does not permit similar static calculations which prompted us to compute the Peierls stress again but now using classical dynamic simulations to facilitate subsequent direct comparisons with RPMD. Another reason for us to perform dynamic simulations was an alarming, even if unnoticed, discrepancy in the literature between “static” and “dynamic” predictions for the Peierls stress reported by different research groups for the same model potential of iron.18,20,25

Shown in Fig. 1 is the flow stress of a single screw dislocation extracted from our MD simulations, plotted as a function of simulation temperature. The flow stress is seen to decrease monotonically with increasing temperature which is typical of dislocations whose motion is temperature- and stress-activated: as the temperature is increased additional thermal energy becomes available helping the dislocation to overcome its motion barriers at an increased rate which requires lower stress to maintain the same dislocation velocity v d . To make sure that our “dynamic” calculations of the Peierls stress are consistent and robust, we performed two series of MD simulations over the same range of temperatures in which the dislocation was forced to move at two different velocities v d : 50 and 10 m/s (see Methods “Molecular dynamics simulations of dislocation motion”). As is seen in Fig. 1, the flow stress at any given temperature is clearly velocity-dependent, however the same Peierls stress is obtained by extrapolation to T = 0 K: τ P = (0.95 ± 0.02) GPa at 50 m/s and τ P = (0.96 ± 0.02) GPa at 10 m/s. Within small error bars, our two “dynamic” predictions agree with our own, as well as previously published “static” predictions for the Peierls stress.18,20

Fig. 1 Classical flow stress of a single screw dislocation in iron predicted in classical MD simulations as a function of temperature. The simulations were performed at two different dislocation velocities v d : 50 m/s (red symbols and red line) and 10 m/s (black symbols and black line). The statistical error bars are smaller than the symbols. The inset depicts the simulation volume in which most atoms are deleted for clarity except for atoms closest to the core of a \({\textstyle{1 \over 2}}\left\langle {111} \right\rangle\) screw dislocation near the center (green) and atoms at the boundaries of the simulation volume (blue and yellow). To persuade the dislocation to move along the horizontal \((11\bar 2)\) plane, the blue atoms in the top and in the bottom layers were moved along the \(\widehat {\bf{x}}\) axis at constant and opposite velocities Full size image

We now turn to RPMD simulations to compute the Peierls stress. Except for the method—classical MD versus quantum RPMD—our classical and quantum simulations were performed using the same interatomic potential and in otherwise identical simulation conditions (see Methods “Ring-polymer molecular dynamics simulations”). Figure 2 presents a comparison between MD and RPMD predictions of flow stress from 5 to 50 K. Following the same extrapolation procedure previously used to extract the Peierls stress from the classical MD data, the Peierls stress predicted in the quantum RPMD simulations is 13% smaller than the classical MD result. In our Supplementary Information 2 we present substantial evidence that our implementation of the RPMD method is highly accurate, see Figs. S2–S4. Thus, we regard our prediction of only a modest reduction in the Peierls stress due to zero-point vibrations in α-iron as more accurate than the dramatic 50% reduction predicted earlier.10

Fig. 2 Temperature dependence of the flow stress of a single screw dislocation predicted in classical (red symbols and red line) and quantum (blue symbols and blue line) simulations. In both cases, the Peierls stress was obtained by a third-order polynomial (solid lines) extrapolation of the finite temperature flow stress to T = 0 K. The orange dashed line is the flow stress as predicted from classical MD simulations at temperature T = 135 K and dislocation velocity v d = 50 m/s Full size image

Although it may be difficult to separate factors that may have contributed to the substantial over-estimation of the effect of zero-point vibrations on the Peierls stress, several explicit or implicit assumptions used in ref.10 to arrive at the dramatic prediction can be suspected. First, reliance on the fully classical TST in defining a critical pathway—a minimum energy path (MEP)—for quantum motion of atoms is inconsistent. Clearly, quantum-mechanical transition states relevant for dislocation motion must be defined and should be searched for in the extended space of closed Feynman paths,26,27 i.e., in 3NP-dimensional RPMD space. Another potential source of error is that Wigner correction (replacing the classical populations of vibrational states with the quantum ones) enhances the population of vibrational modes that accounts for the zero-point energy (ZPE) effects while ignoring other potentially important aspects of quantum motion. That something is amiss about limiting the effect of zero-point vibrations to just the Wigner correction is corroborated by the following qualitative argument.

As shown in Fig. 3, the ZPE of a perfect crystal computed with the same interatomic potential model of iron (see Methods “RPMD calculations of the zero-point energy and the radius of gyration”) is 35 meV/atom. To approximately account for added agitation supplied in the form of zero-point vibrations, let us partition the ZPE equally between the potential and the kinetic energy per atom in the classical system, resulting in a classical temperature of 135 K. By reference to our classical MD results for the flow stress, Fig. 1, this extra temperature would reduce the classical Peierls stress to 0.35 GPa, close to 0.45 GPa reported in ref.10. That our explicit RPMD simulations predict a much more modest reduction suggests that some other effect(s) originating in quantum motion of atoms counterbalances this added agitation. We propose that such a contribution comes from quantum dispersion of atoms, i.e., finite size of a quantum particle compared to a point-particle in classical mechanics.

Fig. 3 Temperature dependence of the per-atom energy in a perfect bcc lattice of α-iron in the quantum (blue symbols and blue line) and the classical (red symbols and red line) systems. Shown in the inset is the excess energy of the quantum system equal to the ZPE at T = 0 K. Shown as a gray vertical band is the Debye temperature of the model of α-iron studied here, T D = (400 ± 20) K Full size image

As was previously observed in the context of quantum diffusion28 and glass transition in quantum liquids,29,30 quantum dispersion can hinder atom rearrangements because in condensed-phase systems the wave-packet of an atom is contracted due to its interaction with neighbor atoms. In other words, each atom is confined/squeezed by its neighboring atoms, which reduces atom dimensions below its free-particle size. In RPMD, how much confinement from its neighbors an atom feels can be assessed by computing the average radius of gyration r G of the ring-polymers representing the atoms. Shown in the inset of Fig. 4 is a plot of r G in a perfect crystal computed for the same model of α-iron as a function of temperature (see Methods “RPMD calculations of the zero-point energy and the radius of gyration”). That atomic confinement is significant can be gauged from another plot in the same figure showing the ratio of r G to the radius of gyration of a free atom \(r_{\mathrm{G}}^{{\mathrm{free}}} = \frac{1}{2} \sqrt {\hbar ^2{\mathrm{/}}mk_{\mathrm{B}}T}\) as a function of temperature. When it comes to thermally activated mechanisms that require overcoming an energy barrier, such as dislocation motion, the already squeezed atoms have to additionally contract their wave packets as they pass through narrow paths leading over the barrier, this additional contraction resulting in an increase in both kinetic and potential energy of atomic configurations corresponding to the barrier states. In the temperature range between 5 and 50 K explored in our RPMD simulations, r G increases by 30% with decreasing temperature partially offsetting the 50% increase in the energy of zero-point vibrations over the same temperature interval (see Fig. 3).

Fig. 4 Temperature dependence of the ensemble-averaged radius of gyration of ring-polymers r G for a perfect bcc lattice of α-iron. \(r_{\mathrm{G}}^{{\mathrm{free}}}\) is the radius of gyration of a free particle of the same mass. The radius of gyration is a measure of the quantum dispersion of the atomic wave function Full size image

As a method for computing equilibrium statistical properties of quantum systems, such as the ZPE and r G , RPMD is asymptotically exact in the limit of continuous Feynman paths (P → ∞). On the other hand, motion dynamics of ring-polymers in RPMD is fictitious and does not necessarily reproduce dynamics of quantum motion of atoms. Every time RPMD is used to study dynamic properties, such as dislocation motion of interest here, the validity of simulation results must be carefully examined. Generally, for an RPMD simulation to be a valid representation of quantum dynamics, the simulated RPMD trajectory should consist of discrete transitions from one potential energy well (basin) to another. Further, the system should reside for a sufficiently long time within each energy basin before transitioning to the next basin, to forget about its preceding dynamic trajectory. This is precisely as expected in the equilibrium quantum TST for which RPMD was recently proven to be exact.26,27,31,32,33,34 For the classical TST to be applicable, the residence time within each basin should exceed the thermalization time. Additionally, when atomic motion is quantum, the same residence time should be longer than the quantum decoherence time (see Supplementary Information 3). In our RPMD simulations, a screw dislocation moves in a stop-and-go fashion at the prescribed average velocity of 50 m/s, as illustrated in Fig. S5 (see Methods “Dislocation motion analysis“). Given the distance between two neighboring energy basins (Peierls valleys) of 2.5 Å in α-iron, this velocity dictates that the dislocation spends on average 5 ps in each Peierls valley, including the residence and the transition segments of its RPMD trajectory as shown in Fig. S6. The quantum thermal time βħ becomes longer than 5 ps below T = 5 K defining the lower bound temperature of our RPMD simulations. The upper bound of 50 K is chosen judiciously, since we observed that at temperatures above 50 K most of the 5 ps time is spent in transition between the Peierls valleys rather than in residence in one of them. To extend the trust range of our RPMD simulations to temperatures above 50 K, the dislocation should be allowed to spend longer time in each Peierls valley necessitating simulations at still lower dislocation velocities.

Our RPMD simulations of quantum motion of dislocations are unprecedented in scale and accuracy. Indeed, all previous RPMD simulations dealt with systems not exceeding a few thousand atoms at temperatures typically above 0.1 of the Debye temperature T D . By comparison, our RPMD simulations entail extended crystal defects embedded in a crystal lattice comprised of ≈150,000 atoms at temperatures as low as 0.01T D . But RPMD method’s rigor comes at a great computational cost: to achieve convergence of quantum RPMD simulations under such conditions, up to 200 replicas of the entire 150,000 atom system had to be integrated using time steps comprising a small fraction of one femtosecond over many nanoseconds. In terms of its raw computational cost, each RPMD simulation takes hundreds and thousands times more compute cycles per unit of simulated time compared to a classical MD simulation with the same interatomic potential. We achieved the needed computational performance by building our RPMD simulations on top of an existing highly efficient MD code, LAMMPS35, and optimizing the resulting RPMD-LAMMPS method for massively parallel computing. To help make the RPMD method more accessible to materials scientists, we make our code freely available and give sufficient technical details of the method’s implementation and fundamentals in our Supplementary Information 1 and 2.

We have developed an efficient implementation of the RPMD method enabling the study of extended defects in materials accounting accurately for quantum corrections to atom dynamics. We employ RPMD for a detailed investigation of the effect of quantum dynamics on dislocation motion in α-iron, which had been previously suggested as a source of the notorious discrepancy between calculated and measured values of the Peierls barrier in bcc metals. The present work utilizes RPMD to show that such corrections are smaller than previously suggested, likely due to quantum dispersion partially compensating the effect of ZPE. The results thus establish that the long-standing discrepancy between calculated and experimental estimates of the Peierls stress in bcc metals remains and further work is needed to settle this outstanding issue in physical metallurgy. We speculate that the origin of this discrepancy is not in the intrinsic properties of single dislocations, rather it lies in complexities emerging from the dislocation network36 that are unavoidable in experimental studies and neglected in most atomistic simulation studies.