Probing the limits of metal plasticity with molecular dynamics simulations

ABSTRACT Ordinarily, the strength and plasticity properties of a metal are defined by dislocations—line defects in the crystal lattice whose motion results in material slippage along lattice

planes1. Dislocation dynamics models are usually used as mesoscale proxies for true atomistic dynamics, which are computationally expensive to perform routinely2. However, atomistic

simulations accurately capture every possible mechanism of material response, resolving every “jiggle and wiggle”3 of atomic motion, whereas dislocation dynamics models do not. Here we

present fully dynamic atomistic simulations of bulk single-crystal plasticity in the body-centred-cubic metal tantalum. Our goal is to quantify the conditions under which the limits of

dislocation-mediated plasticity are reached and to understand what happens to the metal beyond any such limit. In our simulations, the metal is compressed at ultrahigh strain rates along its

[001] crystal axis under conditions of constant pressure, temperature and strain rate. To address the complexity of crystal plasticity processes on the length scales (85–340 nm) and

timescales (1 ns–1μs) that we examine, we use recently developed methods of _in situ_ computational microscopy4,5 to recast the enormous amount of transient trajectory data generated in our

simulations into a form that can be analysed by a human. Our simulations predict that, on reaching certain limiting conditions of strain, dislocations alone can no longer relieve mechanical

loads; instead, another mechanism, known as deformation twinning (the sudden re-orientation of the crystal lattice6), takes over as the dominant mode of dynamic response. Below this limit,

the metal assumes a strain-path-independent steady state of plastic flow in which the flow stress and the dislocation density remain constant as long as the conditions of straining

thereafter remain unchanged. In this distinct state, tantalum flows like a viscous fluid while retaining its crystal lattice and remaining a strong and stiff metal. Access through your

institution Buy or subscribe This is a preview of subscription content, access via your institution ACCESS OPTIONS Access through your institution Access Nature and 54 other Nature Portfolio

journals Get Nature+, our best-value online-access subscription $29.99 / 30 days cancel any time Learn more Subscribe to this journal Receive 51 print issues and online access $199.00 per

year only $3.90 per issue Learn more Buy this article * Purchase on SpringerLink * Instant access to full article PDF Buy now Prices may be subject to local taxes which are calculated during

checkout ADDITIONAL ACCESS OPTIONS: * Log in * Learn about institutional subscriptions * Read our FAQs * Contact customer support SIMILAR CONTENT BEING VIEWED BY OTHERS STATISTICS OF

DISLOCATION AVALANCHES IN FCC AND BCC METALS: DISLOCATION MECHANISMS AND MEAN SWEPT DISTANCES ACROSS MICROSAMPLE SIZES AND TEMPERATURES Article Open access 04 November 2020 IN SITU

ATOMIC-SCALE OBSERVATION OF DISLOCATION CLIMB AND GRAIN BOUNDARY EVOLUTION IN NANOSTRUCTURED METAL Article Open access 18 July 2022 OBSERVING FORMATION AND EVOLUTION OF DISLOCATION CELLS

DURING PLASTIC DEFORMATION Article Open access 13 March 2025 REFERENCES * Bulatov, V. V. & Cai, W. _Computer Simulations of Dislocations_ 196–240 (Oxford Univ. Press, 2006) * Ghoniem, N.

M. et al. Parametric dislocation dynamics: a thermodynamics-based approach to investigations of mesoscopic plastic deformation. _Phys. Rev. B._ 61, 913–927 (2000) Article  ADS  CAS  Google

Scholar  * Feynman, R. _Lectures on Physics_ Vol. 1, 3–6 (1963) ADS  Google Scholar  * Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO — the Open

Visualization Tool. _Model. Simul. Mater. Sci. Eng._ 18, 015012 (2010) Article  ADS  Google Scholar  * Stukowski, A. & Albe, K. Extracting dislocations and non-dislocation crystal

defects from atomistic simulation data. _Model. Simul. Mater. Sci. Eng._ 18, 085001 (2010) Article  ADS  Google Scholar  * Christian, J. W. & Mahajan, S. Deformation twinning. _Prog.

Mater. Sci._ 39, 1–157 (1995) Article  Google Scholar  * Hoge, K. G. & Mukherjee, A. K. The temperature and strain rate dependence of the flow stress of tantalum. _J. Mater. Sci._ 12,

1666–1672 (1977) Article  ADS  CAS  Google Scholar  * Tramontina, D. et al. Molecular dynamics simulations of shock-induced plasticity in tantalum. _High Energy Density Phys._ 10, 9–15

(2014) Article  ADS  CAS  Google Scholar  * Mitchell, T. E. & Spitzig, W. A. Three-stage hardening in tantalum single crystals. _Acta Metall._ 13, 1169–1179 (1965) Article  CAS  Google

Scholar  * Frank, F. C. & Read, W. T. Jr. Multiplication processes for slow moving dislocations. _Phys. Rev._ 79, 722–723 (1950) Article  ADS  CAS  Google Scholar  * Meyers, M. A. et al.

The onset of twinning in metals: a constitutive description. _Acta Mater._ 49, 4025–4039 (2001) Article  CAS  Google Scholar  * Shields, J. A. et al. Deformation of high purity tantalum

single crystals at 4.2 K. _Mater. Sci. Eng._ 20, 71–81 (1975) Article  CAS  Google Scholar  * Cotterill, R. M. J. Does dislocation density have a natural limit? _Phys. Lett. A_ 60, 61–62

(1977) Article  ADS  Google Scholar  * Florando, J. N. et al. Effect of strain rate and dislocation density on the twinning behavior in tantalum. _AIP Adv._ 6, 045120 (2016) Article  ADS 

Google Scholar  * Sleeswyk, A. W. 1/2〈111〉 screw dislocations and nucleation of {112}〈111〉 twins in the b.c.c. lattice. _Phil. Mag._ 8, 1467–1486 (1963) Article  ADS  Google Scholar  *

Marian, J., Cai, W. & Bulatov, V. V. Dynamic transitions from smooth to rough to twinning in dislocation motion. _Nat. Mater._ 3, 158–163 (2004) Article  ADS  CAS  Google Scholar  *

Martyushev, L. M. & Seleznev, V. D. Maximum entropy production principle in physics, chemistry and biology. _Phys. Rep._ 426, 1–45 (2006) Article  ADS  MathSciNet  CAS  Google Scholar  *

Hsiung, L. L. Shock-induced phase transformation in tantalum. _J. Phys. Condens. Matter_ 22, 385702 (2010) Article  ADS  Google Scholar  * Nemat-Nasser, S. et al. Microstructure of

high-strain, high-strain-rate deformed tantalum. _Acta Mater._ 46, 1307–1325 (1998) Article  CAS  Google Scholar  * Lu, C. H. et al. Phase transformation in tantalum under extreme laser

deformation. _Sci. Rep._ 5, 15064 (2015) Article  ADS  CAS  Google Scholar  * Reed, B. W. et al. A unified approach for extracting strength information from non-simple compression waves.

Part II. Experiment and comparison with simulation. _J. Appl. Phys._ 110, 113506 (2011) Article  ADS  Google Scholar  * Johnston, W. G. & Gilman, J. J. Dislocation velocities,

dislocation densities and plastic flow in lithium fluoride crystals. _J. Appl. Phys._ 30, 129–144 (1959) Article  ADS  CAS  Google Scholar  * Mecking, H. & Kocks, U. F. Kinetics of flow

and strain hardening. _Acta Metall._ 29, 1865–1875 (1981) Article  CAS  Google Scholar  * Seeger, A. Evidence of enhanced self-organization in the work-hardening stage V of fcc metals.

_Philos. Mag. Lett._ 81, 129–136 (2001) Article  ADS  CAS  Google Scholar  * Banerjee, J. K. Barreling of solid cylinders under axial compression. _J. Eng. Mater. Technol._ 107, 138–144

(1985) Article  Google Scholar  * Roylance, D. _Stress–Strain Curves_

https://ocw.mit.edu/courses/materials-science-and-engineering/3-11-mechanics-of-materials-fall-1999/modules/MIT3_11F99_ss.pdf (MIT OpenCourseWare, 2001) * Khan, A. S. & Huang, S.

_Continuum Theory of Plasticity_ 37–40 (Wiley-Interscience, 1995) * Carpay, F. M. A. et al. Constrained deformation of molybdenum single crystals. _Acta Metall._ 23, 1473–1478 (1975) Article

  CAS  Google Scholar  * Saada, G. On hardening due to the recombination of dislocations. _Acta Metall._ 8, 841–847 (1960) Article  CAS  Google Scholar  * Bulatov, V. V. et al. Dislocation

multi-junctions and strain hardening. _Nature_ 440, 1174–1178 (2006) Article  ADS  CAS  Google Scholar  * Louchet, F. & Viguier, B. Ordinary dislocations in γ-TiAl: cusp unzipping, jog

dragging and stress anomaly. _Philos. Mag. A_ 80, 765–779 (2000) Article  ADS  CAS  Google Scholar  * Christian, J. W. & Masters, B. C. Low-temperature deformation of body-centered cubic

metals. _Proc. R. Soc. Lond. A_ 281, 223–239 (1964) Article  ADS  CAS  Google Scholar  * Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. _J. Comput. Phys._ 117,

1–19 (1995) Article  ADS  CAS  Google Scholar  * Li, Y. H. et al. Embedded-atom-method tantalum potential developed by force-matching method. _Phys. Rev. B_ 67, 125101 (2003) Article  ADS 

Google Scholar  * Boxall, A. IBM’s Sequoia tops the world’s fastest supercomputer list. _Digital

Trends_http://www.digitaltrends.com/computing/ibms-sequoia-tops-the-worlds-fastest-supercomputer-list/ (2012) * Nguyen, L. D. & Warner, D. H. Improbability of void growth in aluminum via

dislocation nucleation under typical laboratory conditions. _Phys. Rev. Lett._ 108, 035501 (2012) Article  ADS  CAS  Google Scholar  * Mauro, J. C. & Smedskjaer, M. M. Unified physics

of stretched exponential relaxation and Weibull fracture statistics. _Physica A_ 391, 6121–6127 (2012) Article  ADS  Google Scholar  * Hähner, P. et al. Fractal dislocation patterning during

plastic deformation. _Phys. Rev. Lett._ 81, 2470–2473 (1998) Article  ADS  Google Scholar  Download references ACKNOWLEDGEMENTS This work was performed under the auspices of the US

Department of Energy by Lawrence Livermore National Laboratory under contract W-7405-Eng-48. This work was supported by the NNSA ASC programme. Computing support for this work came from the

Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge programme and Jülich Supercomputing Centre at Forschungszentrum Jülich, Germany. AUTHOR INFORMATION

AUTHORS AND AFFILIATIONS * Lawrence Livermore National Laboratory, Livermore, California, USA Luis A. Zepeda-Ruiz, Tomas Oppelstrup & Vasily V. Bulatov * Technische Universität

Darmstadt, Darmstadt, Germany Alexander Stukowski Authors * Luis A. Zepeda-Ruiz View author publications You can also search for this author inPubMed Google Scholar * Alexander Stukowski

View author publications You can also search for this author inPubMed Google Scholar * Tomas Oppelstrup View author publications You can also search for this author inPubMed Google Scholar *

Vasily V. Bulatov View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS L.A.Z.-R. ran most of the molecular dynamics simulations and analysed

the results, A.S. ran molecular dynamics simulations and developed methods for computational microscopy and visualization, T.O. optimized run-time efficiency and data management of molecular

dynamics simulations, and V.V.B. developed the concept, planned the research, generated starting configurations for molecular dynamics simulations, analysed the results and wrote the paper.

CORRESPONDING AUTHOR Correspondence to Vasily V. Bulatov. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing financial interests. ADDITIONAL INFORMATION REVIEWER

INFORMATION _Nature_ thanks M. Zaiser and the other anonymous reviewer(s) for their contribution to the peer review of this work. Publisher's note: Springer Nature remains neutral with

regard to jurisdictional claims in published maps and institutional affiliations. EXTENDED DATA FIGURES AND TABLES EXTENDED DATA FIGURE 1 PLASTIC YIELD RESPONSE DEPENDS ON THE INITIAL

DENSITY OF DISLOCATION SOURCES. Stress as a function of true strain and specimen size computed in three MD simulations of compression at rate ×25 from three different initial configurations

of dislocation sources. EXTENDED DATA FIGURE 2 DETECTING TWINNING DURING STRAINING SIMULATIONS. Stress (dashed curve) and the volume fraction of twins (solid curve) as a function of true

strain under compression at rate ×5 and temperature 25 K. The twin fraction was computed using the GSA. EXTENDED DATA FIGURE 3 CONTINUOUS COOLING DURING STRAINING PERMITS DETECTION OF A

TWINNING TRANSITION. Volume fraction of twins (solid curve) and flow stress (dashed curve) as a function of temperature, computed in a simulation at fixed rate ×5 in which temperature was

reduced at a constant rate from 300 K to 10 K. A twinning transition is identified by the temperature at which the twin fraction begins to rise rapidly from zero. Preceding this simulation,

the crystal was pre-strained at the same rate ×5 and a fixed temperature _T_ = 300 K, where it attained a steady flow stress of 3.2 GPa. EXTENDED DATA FIGURE 4 RELAXATION OF DISLOCATION

DENSITY AFTER UNLOADING OF TWO CRYSTALS PRE-STRAINED AT RATE ×75. The blue line represents isothermal relaxation after isothermal straining at _T_ = 300 K. The red line depicts additional

relaxation after two opposite surfaces of the simulated crystal were exposed to vacuum. The black line is adiabatic relaxation after adiabatic straining. Both isothermal and adiabatic

relaxation simulations start at 2.5 ns in the end of isothermal and adiabatic pre-straining simulations (not shown). EXTENDED DATA FIGURE 5 KNEADING THE METAL. Stress as a function of true

(von Mises) strain, computed under compression at constant true rate ×25 along the three principal axes of the crystal. After compressing the crystal to one-quarter of its initial length

along the _z_ axis, the strain axis is changed from _z_ to _y_, from _y_ to _x_, and then from _x_ back to _z_. Letters above the stress–strain curves label the axes for each compression

cycle. EXTENDED DATA FIGURE 6 EVOLUTION OF DISLOCATION NETWORK TOPOLOGY UNDER COMPRESSIVE STRAIN AT RATE ×1. Following rapid dislocation multiplication at yield, regular binary junctions

appear first (red line) closely followed by ternary multi-junctions (black line). After reaching stationary flow at a strain of about 0.4, dislocation density and network composition

(topology) remain stationary within statistical noise. EXTENDED DATA FIGURE 7 NETWORK EVOLUTION ALONG A STEPWISE INCREASE IN THE STRAIN RATES. The number of binary junctions along a stepwise

strain trajectory is shown as red circles and the solid red line, and the number of ternary junctions is shown as green squares and the solid green line. The dashed blue line shows the

stepwise strain trajectory with strain rates marked at each rate step. The thin solid lines show the numbers of binary (red) and ternary (green) junctions along continuations of the

interrupted straining steps. The inset shows the ratio of the number of ternary junctions to the number of binary junctions attained in the saturated flow state as a function of strain rate.

Error bars are the standard deviation from the mean values. EXTENDED DATA FIGURE 8 EVOLUTION OF DISLOCATION CHARACTERS UNDER COMPRESSION AT RATE ×25. A, Histograms of dislocation character

angle distributions computed for configurations A (black bar near zero), B (red bars) and C (open bars), which are marked along the strain trajectory in B. For reference, the dashed line

depicts a hypothetical uniform distribution of character angles with the same integral length of dislocations as in configuration C. The histogram counts are over the bins along the

log[cos(_θ_)] axis. B, Ratio of the total length of near-screw dislocations to the total length of near-edge dislocations as a function of strain. EXTENDED DATA FIGURE 9 DIFFERENTIAL SLIP

TRACE ANALYSIS REVEALS HOW DISLOCATIONS MOVE AND INTERACT. A, The blue lines show positions and shapes of the dislocation lines in the initial configuration and the green lines show

dislocations in the final configuration attained a few picoseconds later. The grey ‘slip traces’ consist of atoms whose local von Mises shear strain accumulated between the initial and the

final dislocation positions exceeds the 0.15 threshold. Only a relatively small fraction of dislocations had swept substantial areas, whereas positions of most other dislocations in the two

configurations coincide, suggesting little or no motion over the time interval. B, A magnified fragment of the same differential plot showing grey areas swept by several dislocations in more

detail. The smooth step seen on the slip trace in the foreground reveals a ‘jog’ (a turn of a dislocation line inside a crystal) on the moving dislocation. C, Cross-slip of a screw

dislocation from its initial position (blue) to its final position (green). The shape of the cylindrical traced area reveals the detailed trajectory of the screw dislocation between its

initial and final positions. D, Annihilation of two dislocations, as evidenced by a slip trace area bounded on its two sides by two blue lines: one straight screw dislocation above and one

curved dislocation below. That annihilation has taken place is deduced from the absence of green lines, which would otherwise show final positions of the two dislocations. EXTENDED DATA

FIGURE 10 EFFECT OF SIMULATION VOLUME SIZE ON THE STRAIN RESPONSE. A, Stress–strain response under straining at rate ×1 (1.11 × 107 s−1). The thick red line is the stress–strain response

simulated in a volume eight times greater (about 268 million atoms) than the one used in most other simulations. The thin grey lines correspond to eight independent simulations at the same

rate, but performed in the standard-sized volume with about 33 million atoms. The thick black line was obtained by averaging over these eight simulations. B, The corresponding density–strain

curves, with line colours and types matching the stress–strain curves on the left. SUPPLEMENTARY INFORMATION DISLOCATION MULTIPLICATION FROM THE INITIAL SOURCES RESULTS IN THE DEVELOPMENT

OF A DENSE DISLOCATION NETWORK. Crystal containing dislocations sources (loops) is subjected to uniaxial compression along the [001] axis at a constant true straining rate of 2.78.108 1/s

(x25). The simulation volume contains about 268 million atoms of tantalum. The video sequence progresses through extension of the initial hexagon-shaped loops, to dislocation collisions

resulting in the formation of dislocation junctions, to an increasingly dense dislocation network. Dislocation positions, shapes and Burgers vectors were extracted using the DXA algorithm5.

All atoms and defects other than dislocations, such as vacancies, interstitials and clusters, are removed for clarity. The green lines represent ½<111> dislocations and the pink lines

depict <100> junction dislocations. (MP4 27123 kb) CRYSTAL MICROSTRUCTURE EVOLUTION UNDER STRAINING AT RATE X50. In this MD simulation a crystal containing dislocations sources (loops)

was subjected to uniaxial compression along the [001] axis at a constant true straining rate of 5.56.108 1/s (x50). The simulation volume contains about 33 million atoms of tantalum. This

video sequence progresses through extension of the initial loops, to nucleation of embryonic twins on screw dislocations, to rapid propagation and growth of twinning particles. The meaning

of colours is as defined in the caption to Fig. 1: the outer surfaces bounding the twins are coloured light grey whereas the insides of twin particles are coloured red, yellow, magenta or

cyan depending to each twin's rotational variant. (MP4 19343 kb) “METAL KNEADING” AT RATE X25. This MD simulation was performed on a brick-shaped tantalum crystal with the ratio of

initial box dimensions 1:2:4. After full compression along Z axis to ¼ of its initial dimension the brick’s shape becomes 2:4:1 (due to Poisson’s expansion in two lateral dimensions, the

brick’s volume remains very nearly constant under compression) another MD simulation starts in which the brick is compressed along its now longest Y-axis. After the second compression cycle

is completed, the brick is compressed along its now longest X-axis. After three compression cycles the brick recovers its initial shape 1:2:4 and one more Z-axis compression cycle is

performed (see related stress-strain plots in Extended Data Figure 5). (MP4 26846 kb) DISLOCATION MOTION IN MORE DETAIL This simulation was performed at rate 1.11.107 1/s (x1) from a

configuration attained past yield under pre-straining at rate 5.55.107 1/s (x5). Reduction in dislocation density can be observed over the first few frames immediately following the sudden

drop in the straining rate from x5 to x1 at time _t_=0. Subsequently the network attains a dynamic steady state in which dislocation multiplication is balanced by dislocation annihilation.

Taken at more frequent time intervals, this sequence reveals various events in the life of dislocations in greater detail than in the other videos. One can observe that dislocation motion is

not steady but proceeds in a stop-and-go manner which is also revealed in Extended Data Fig. 9a. (MP4 26951 kb) POWERPOINT SLIDES POWERPOINT SLIDE FOR FIG. 1 POWERPOINT SLIDE FOR FIG. 2

POWERPOINT SLIDE FOR FIG. 3 POWERPOINT SLIDE FOR FIG. 4 POWERPOINT SLIDE FOR FIG. 5 RIGHTS AND PERMISSIONS Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Zepeda-Ruiz, L.,

Stukowski, A., Oppelstrup, T. _et al._ Probing the limits of metal plasticity with molecular dynamics simulations. _Nature_ 550, 492–495 (2017). https://doi.org/10.1038/nature23472 Download

citation * Received: 28 April 2016 * Accepted: 23 June 2017 * Published: 27 September 2017 * Issue Date: 26 October 2017 * DOI: https://doi.org/10.1038/nature23472 SHARE THIS ARTICLE Anyone

you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the

Springer Nature SharedIt content-sharing initiative