Quantum effects on dislocation motion from Ring-Polymer Molecular Dynamics

Quantum motion of atoms known as zero-point vibrations is recognized to be important at low temperatures in condensed matter systems comprised of light atoms or ions, affecting such properties and behaviors as proton-transfer reactions, vibrational spectra of water and ice, and mechanical properties of low temperature helium. Recently, quantum motion of atoms was proposed to explain a long-standing discrepancy between theoretically computed and experimentally measured low-temperature resistance (Peierls stress) to dislocation motion in iron and possibly other metals with high atomic masses. Here we report the first direct simulations of quantum motion of screw dislocations in iron within the exact formalism of Ring-Polymer Molecular Dynamics (RPMD) that rigorously accounts for quantum effects on the statistics of condensed-phase systems. Our quantum RPMD simulations predict only a modest ($\approx\!13\%$) reduction in the Peierls stress in iron compared to its fully classical prediction. Our simulations confirm that reduction in the Peierls stress solely due to the zero-point energy is close to $50\%$ predicted earlier, but its effect is substantially offset by an increase in the effective atom size with decreasing temperature, an effect known as quantum dispersion. Thus, quantum motion of atoms does not resolve the notorious discrepancy between theoretical and experimental values of the Peierls stress in iron.


Introduction
Quantum motion of atoms known as zero-point vibrations is recognized to be important at low temperatures in condensed matter systems comprised of light atoms or ions, affecting such properties and behaviors as protontransfer reactions 1,2 , vibrational spectra of water 3,4 and ice 5,6 , and mechanical properties of low temperature helium [7][8][9] . Recently, quantum motion of atoms was proposed to explain a long-standing discrepancy between theoretically computed and experimentally measured low-temperature resistance (Peierls stress) to dislocation motion in iron and possibly other metals with high atomic masses [10][11][12] .
Experimental estimates for the Peierls stress τ P in body-centered cubic (bcc) metals are obtained by relating this parameter to the low-temperature yield strength of the metal extrapolated to zero temperature. In α-iron, the so-obtained experimental estimates for τ P range between 0.35 GPa and 0.45 GPa [13][14][15][16] whereas atomistic calculations consistently predict τ P to be much higher, ranging between 0.9 GPa and 1.2 GPa depending on the interatomic potential model used in the calculations 17,18,10 . At first perceived as a failure of interatomic potentials, subsequent electronic structure calculations confirmed the stubbornly high theoretical values of τ P in α-iron to be between 1.0 GPa and 1.4 GPa [19][20][21] . In ref. 10 the effect of zero-point vibrations was accounted for by computing an approximate quantum correction to the rate of dislocation motion predicted within the classical transition state theory (TST) 10 . The resulting correction reduced the Peierls stress from its fully classical value of 0.9 GPa by about 50%, to τ P = 0.45 GPa thus bringing it in close agreement with the experimental estimates.
Although appealing, the dramatic two-fold reduction in the Peierls stress was predicted from a harmonic approximation to the classical TST in which classical populations of vibrational modes were replaced with their quantum analogs. Given its potential importance, here we re-examine the effect of zero-point vibrations on dislocation motion in α-iron using Ring-Polymer Molecular Dynamics 22,23 simulations (RPMD). Based on Feynman's path integral formulation of quantum mechanics, RPMD is a method that fully accounts for all effects of quantum motion of atoms on equilibrium statistical properties. As detailed in fig. S1, RPMD amounts to a simultaneous simulation of P replicas of the entire N-atom system in which each atom is represented by its P clones connected into a "polymer ring" (a closed Feynman path) by elastic springs. RPMD is asymptotically exact in the limit P → ∞. To take full advantage of RPMD accuracy, in the following we make no assumptions on the character of dislocation motion and measure the effect of zero-point vibrations on the Peierls stress in the most direct manner possible, akin to the very method by which the Peierls stress is estimated in experiment. Using the same interatomic model of α-iron as in ref. 10, here we compute the resistance to dislocation motion twice -first using fully classical Molecular Dynamics (MD) and then using RPMD simulations -and compare the results in the limit of zero temperature to assess the differences.

Results and Discussion
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 1 2 111 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 potential 10,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 same 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 m/s and 10 m/s (see Methods 1). 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 . 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 2). Figure 2 presents a comparison between MD and RPMD predictions of flow stress from 5 K 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, 3 S3, and 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 . 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 3) 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.  As was previously observed in the context of quantum diffusion 28 and glass transition in quantum liquids 29,30 , quantum dispersion can hinder atom re-arrangements 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 sees 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 3). 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 free G = 2 /mk B T /2 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 K 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). 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 4). 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.01 T 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 LAMMPS 35 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 network 36 that are unavoidable in experimental studies and neglected in most atomistic simulation studies.

Molecular Dynamics simulations of dislocation motion
All atomistic simulations were performed in a volume shaped as a rectangular prism oriented as shown in the inset in fig. 1  Classical MD simulations were performed using LAMMPS 35 code with an embedded-atom method (EAM) potential previously developed for the bcc α-phase of elemental iron 39 . The simulation timestep was ∆t = 3 fs, ≈ 1/40 th of the period of the highest frequency mode in the phonon spectrum computed for this model system. The Debye temperature, T D = (400 ± 20) K, was estimated from the same phonon spectrum, as described in ref. 40. Temperature was maintained near-constant (NVT ensemble) in each simulation using a Langevin thermostat with the temperature relaxation time of 30 ps, which was found to be sufficiently gentle not to affect the flow stress when compared to corresponding microcanonical simulations. Classical NVT simulations were performed at temperatures ranging from 0.1 K to 300 K using per-atom volumes adjusted at each temperature to maintain internal pressure close to zero. Prior to running mobility simulations, the model was thermalized at the desired temperatures for 300 ps after which two rigid layers (slabs) of atoms at the top and the bottom of the simulation volume (blue atoms in fig. 1) were displaced at a constant velocity in the opposite directions along thex-axis ([111]). The slab velocity was adjusted so as to make the dislocation move at a certain average velocity (50 m/s or 10 m/s) once its motion becomes steady after an initial transient period. The straining rate was kept constant for 10 ns over which the τ xz component of the internal stress was recorded as a measure of the flow stress. A stress versus time curve obtained in one such MD simulation is shown in fig. S7.

Ring-Polymer Molecular Dynamics simulations
The RPMD simulations were performed using our own implementation of the method in LAMMPS (Supplementary Information 1 and 2). Except for the significantly greater number of degrees of freedom -3N in MD and 3N P in RPMD -and a shorter time step, all other details of RPMD simulations were the same as in the classical MD simulations described in the preceding section (Methods 1).
The number of beads P (system's replicas) was selected individually for each temperature based on convergence of the system's energy and internal stress-tensor with the number of beads. Figure S8 presents data obtained in one such convergence study at T = 30 K. Based on this data we elected to use 40 beads for RPMD simulations at this particular temperature. Based on similar analyses we ended up using 200, 100, 50, 40, 30, and 20 beads at temperatures 5 K, 10 K, 20 K, 30 K, 40 K, and 50 K, respectively. The timestep was set equal to 1/25th of the period of the highest frequency eigenmode of the ring-polymers 41 , which in practice resulted in time steps of 0.48 fs for 5 K, 10 K, 20 K, and 30 K, and 0.40 fs for 40 K and 50 K. Similar to the classical MD, we adjusted the per-atom volume at each temperature in the quantum system to maintain its internal pressure close to zero; in practice the difference in the zero-pressure lattice parameter between MD and RPMD was only of about ±0.05% (at the same temperature) and had no appreciable effect on dislocation motion.
To induce dislocation motion and to compute the Peierls stress we used the same procedures as in the classical MD simulations. To additionally verify that our RPMD predictions are accurate, we repeated RPMD simulations at T = 40 K and 50 K using twice as many beads as was judged necessary for fully converged simulations and observed no appreciable changes in the results. Additionally, we verified that our flow stress predictions are converged with respect to time step by repeating our RPMD simulations at T = 30 K with time steps of 0.20 and 0.10 fs, i.e. 1 2 and 1 4 of the 0.40 fs deemed acceptable at this temperature.

RPMD calculations of the zero-point energy and the radius of gyration
RPMD simulations of a perfect bcc lattice of iron with 5, 488 atoms were used to compute the ZPE and r G shown in figs. 3 and 4. Simulations were performed at the same temperatures and with the same numbers of beads P as in Methods 2. Additionally, RPMD simulations were performed at higher temperatures: 100 K (with 16 beads), 150 K (12 beads), 200 K (12 beads), 250 K (10 beads), 300 K (9 beads), 400 K (8 beads), 500 K (7 beads), and 600 K (7 beads). The classical energies shown on the plot were extracted from MD simulations of the same system at the same temperatures. The radius of gyration r G of the ring-polymers in RPMD shown in fig. 4 was computed as where . . . denotes the canonical ensemble average and is the ring-polymer centroid. The radius of gyration of the free atom shown in the same figure was computed analytically 42 :

Dislocation motion analysis
For RPMD simulations to be representative of stochastic dynamics of dislocation motion, it was necessary to assess whether the dislocation spends sufficient time in its energy basins between basin-to-basin transitions. The analysis was performed on the snapshots of atom centroid positions, eq. (2), saved every 1.4 fs along the simulated RPMD trajectories. Within each time snapshot the dislocation's position was taken as a center of mass of atom centroids that were deemed "defective" according to the adaptive Common Neighbor Analysis 43,44 ; in our low-temperature RPMD simulations most such "defective" atoms belong to the dislocation core. Figure S6a is an example of a dislocation trajectory extracted from one of our RPMD simulations. The trajectory reveals a stop-and-go character of dislocation motion in which periods of residence in the energy basins (horizontal segments) are interrupted by more or less fast transitions between the basins (step segments). Duration of the residence periods decreases with increasing temperature ( fig. S6b) and at temperatures above 50 K distinct residence periods are no longer recognized. Short of running RPMD simulations at a lower dislocation velocity (which would extend the average residence periods but also raise the computational cost) we no longer consider our RPMD simulations at temperatures above 50 K to be representative of dislocation motion in our model of α-iron. 9