Machine-learning molecular dynamics on the LUMI supercomputer reaches the 2-MeV regime with billion-atom simulations, revealing how tungsten responds to the extreme particle bombardment expected inside future fusion reactors.
For decades, scientists have understood the basic mechanism by which energetic particles damage metals: a high-energy particle strikes an atom, knocking it from its lattice site and triggering a rapidly expanding collision cascade. But understanding the first few trillionths of a second of that event in a material as important as tungsten is considerably harder than the basic description suggests.
Now, researchers from the University of Helsinki, Åbo Akademi University and CSC–IT Center for Science have used machine-learning-driven molecular dynamics at unprecedented scale to follow radiation-damage cascades in tungsten from just 40 electronvolts to 2 megaelectronvolts.
The calculations reached systems containing up to one billion tungsten atoms and were performed on the GPU nodes of Europe's LUMI supercomputer. The simulations reveal, for the first time, a complete progression through four distinct regimes of primary radiation damage, including a previously inaccessible high-energy regime in which defect production returns to a linear relationship with deposited energy.
It is a remarkable demonstration of what happens when machine learning, GPU computing, and molecular dynamics are combined at extreme scale.
And for fusion research, the result could provide a more accurate computational foundation for predicting how reactor materials deteriorate under neutron bombardment.
Tungsten meets the fusion environment
Tungsten is one of the leading candidates for the plasma-facing components of future fusion reactors.
Its appeal is straightforward: it has an exceptionally high melting point and can withstand extreme thermal and radiation environments.
But the same fusion reactions that produce energy also create an extraordinarily hostile particle environment.
A deuterium-tritium fusion reaction produces a 14.1-MeV neutron. When those neutrons strike tungsten, they can transfer as much as approximately 300 keV of recoil energy to tungsten atoms.
That recoil initiates a cascade.
One displaced atom strikes another.
That atom displaces another.
Within an incredibly short period, hundreds, thousands, or potentially many more atoms can be pushed away from their normal lattice positions.
The resulting vacancies, self-interstitial atoms, dislocation loops, and defect clusters ultimately contribute to swelling, embrittlement, and other forms of material degradation.
For fusion engineers, predicting primary damage accurately is essential.
For computational scientists, however, there is a problem.
The cascade is enormous.
The computational wall
Molecular dynamics is one of the most powerful tools available for studying radiation damage because it follows atoms individually according to the underlying interatomic forces.
But that precision comes at a price.
At relatively low recoil energies, a simulation containing thousands or millions of atoms can be sufficient.
At higher energies, the collision cascade becomes physically larger.
A simulation box that is too small causes the cascade to interact with its own boundaries, contaminating the physics.
Before this work, full atomistic simulations had reached approximately 300 keV, enough to cover the fusion-neutron energy range but leaving the behavior at MeV energies largely unexplored computationally.
Experiments, meanwhile, routinely investigate MeV-scale recoil energies using heavy-ion irradiation.
That created a significant gap between what experimentalists could produce and what computational scientists could model atom by atom.
The new work attacks that gap directly.
One billion atoms at a time
The researchers simulated primary knock-on atom, or PKA, energies spanning six orders of magnitude, from 40 eV to 2 MeV.
At the highest energies, the simulation cells contained approximately 1.024 billion tungsten atoms in a cube about 255 nanometers on each side. Twenty independent simulations were performed at each of the 1- and 2-MeV energies.
The progression in system size is striking.
At 50 keV, the calculations used about 8.2 million atoms.
At 100 keV, approximately 16 million.
At 200 keV, 54 million.
At 300–500 keV, 128 million.
And at 1–2 MeV:
more than one billion atoms.
This is not simply a bigger simulation.
It represents a fundamentally different computational regime.
Machine learning makes the scale possible
The breakthrough depends heavily on the interatomic potential.
Traditional molecular-dynamics calculations require an accurate description of the forces between atoms while evaluating those forces millions or billions of times during a simulation.
The researchers used tabGAP, a tabulated Gaussian Approximation Potential, a machine-learned interatomic potential designed to provide high computational efficiency while retaining the accuracy needed for atomistic materials modeling.
They ported the LAMMPS implementation of tabGAP to GPUs using the Kokkos performance-portability framework.
This is where the study becomes particularly interesting from an HPC perspective.
The researchers weren't simply given a larger computer.
They redesigned the computational workload to exploit modern accelerator hardware.
The force calculation was parallelized over atoms, and the GPU implementation achieved a 1.6× speedup on a single AMD MI250X GPU GCD compared with the original CPU implementation running on a full LUMI CPU node containing 128 AMD EPYC 7763 cores.
The electronic-stopping calculation used during the high-energy cascades was also ported to Kokkos, allowing the frequently executed cascade calculations to run on GPUs.
This is a powerful example of modern HPC optimization:
better physics + better algorithms + accelerator computing = a previously inaccessible simulation regime.
Following a cascade for 50 picoseconds
The simulations begin with a tungsten lattice relaxed to 300 kelvin and zero pressure.
A primary recoil is then launched into the material.
Because the collision evolves extremely rapidly, the researchers use an adaptive timestep and follow the cascade for 50 picoseconds.
A thin, 8-ångström boundary region is thermostatted at 300 K using a Nosé-Hoover thermostat. This boundary treatment removes heat and damps pressure waves so that the cascade does not artificially interact with the simulation boundary.
At high recoil energies, another physical effect becomes important: the energetic atoms can lose energy to electronic excitations.
The researchers model that electronic stopping as a friction force for atoms above 10 eV using stopping data from SRIM.
Every detail matters.
When a billion-atom calculation is being used to make predictions about a fusion reactor, numerical artifacts can be as dangerous as missing physics.
Turning atoms into data
The raw molecular-dynamics trajectories contain enormous quantities of information.
The researchers needed to determine which atoms had been displaced, where vacancies and self-interstitial atoms formed, how defects clustered and whether dislocations developed.
They used several computational analysis techniques:
* Wigner-Seitz analysis to identify vacancies and self-interstitial atoms;
* the Dislocation Extraction Algorithm (DXA) to identify dislocations;
* cluster analysis to determine defect-cluster sizes; and
* OVITO for analysis and visualization.
The resulting dataset allows the researchers to move beyond simply asking how many atoms were displaced.
They can investigate the morphology of the cascade.
That turns out to be crucial.
Four regimes hidden inside the cascade
The simulations reveal four distinct regimes of radiation damage.
Regime I: Near the displacement threshold
At very low energies, the recoil may barely have enough energy to permanently displace atoms.
The minimum tungsten displacement energy is approximately 42 ± 1 eV for certain crystallographic directions, while the average threshold over directions is approximately 95 eV and the maximum exceeds 250 eV.
Interestingly, the simulations show that around 80–130 eV, only about 0.25 Frenkel pairs are produced per recoil on average.
Even when an atom is displaced, many defects subsequently recombine.
Regime II: The sublinear heat-spike regime
As recoil energy increases, conventional collision cascades generate increasingly intense local heating.
The cascade forms a microscopic heat spike, a dense, transient region that behaves somewhat like a tiny volume of hot liquid.
This promotes recombination.
Consequently, the number of surviving defects grows more slowly than the deposited energy.
The simulations agree well with established arc-dpa-based modeling in this regime.
But then something unexpected happens.
Regime III: The superlinear regime
At approximately 20–30 keV, the behavior changes.
Subcascades begin to form, but they don't necessarily separate cleanly.
Instead, many remain close enough for their heat spikes to overlap.
The result can be an extraordinarily dense region of energy deposition.
The simulations show that these compact or overlapping cascades create anomalously large defect clusters and more surviving defects than predicted by the conventional models.
The researchers identify a heat-spike radius of roughly 3 nanometers at the transition.
This produces a superlinear increase in damage.
In other words, adding more energy doesn't simply produce proportionally more damage.
Under these conditions, the cascade becomes unusually efficient at producing persistent defects.
The computer recreates a microscopic explosion
The visualization of these simulations is extraordinary.
At 20 keV, the cascade remains relatively compact.
At 200 keV, overlapping subcascades become apparent.
At 2 MeV, the simulation reveals a much larger structure in which multiple subcascades evolve and eventually separate.
The researchers analyze the cascades at femtosecond and picosecond timescales, capturing both the initial energetic collisions and the later heat-spike evolution.
The computational scale is enormous, but the physical event itself is fleeting.
The entire primary-damage process unfolds in a fraction of a nanosecond.
Supercomputing effectively provides a microscope for time as well as space.
Regime IV: The Linear Frontier
The biggest discovery comes above approximately 300 keV.
At these energies, subcascades increasingly separate far enough that they no longer overlap.
The dense heat spikes responsible for the superlinear regime stop becoming progressively more extreme.
Instead, the energy is divided among increasingly independent subcascades.
The defect production therefore returns to a linear trend.
This is the first time the high-energy linear regime has been directly revealed and quantified in atomistic tungsten simulations.
And the energy is significant.
Approximately 300 keV is also the maximum recoil energy tungsten can receive from a 14.1-MeV fusion neutron.
That coincidence is extremely useful for fusion research.
It means the new simulations identify the transition right at the upper edge of the primary radiation-damage regime most directly relevant to a fusion reactor.
Why the transition matters
The distinction between the four regimes isn't merely academic.
Most engineering models need to convert radiation energy into an estimate of the number of defects created.
If the damage is assumed to increase linearly when it is actually superlinear, defect production could be underestimated.
If the superlinear behavior is incorrectly extrapolated indefinitely, it could instead be overestimated at higher energies.
The new simulations show that neither assumption is correct across the full energy range.
There is a superlinear window.
Then, as the subcascades separate, the physics changes again.
A new full-range damage model
The researchers use the computational results to construct a revised analytical model covering all four regimes.
The model combines the established arc-dpa formulation with a new energy-dependent enhancement function.
That enhancement function rises through the superlinear regime and then saturates as the system approaches the high-energy linear regime.
The resulting model reproduces the simulated trend from near-threshold energies through the sublinear and superlinear regimes and into the newly observed linear high-energy regime.
This is where an extreme-scale simulation becomes useful to engineers.
The goal isn't to run a billion-atom calculation every time someone wants to estimate radiation damage in a reactor component.
The goal is to use those simulations to build better reduced-order models that can be incorporated into larger materials and reactor simulations.
Billion atoms, but only 50 picoseconds
One of the most fascinating aspects of the research is the mismatch between spatial and temporal scale.
The highest-energy simulations contain one billion atoms.
Yet each cascade is followed for only about 50 picoseconds.
That's because primary radiation damage happens extraordinarily quickly.
The simulation therefore represents a massive three-dimensional computational domain evolving over an almost unimaginably short interval.
This is exactly the kind of workload that modern supercomputers are uniquely suited to handle.
Thousands or millions of atom interactions must be calculated repeatedly while maintaining enough spatial resolution to prevent the cascade from interacting artificially with the boundaries.
Artificial heat spikes confirm the physics
The team also performed controlled simulations in which kinetic energy was artificially deposited into a spherical region of tungsten.
These experiments effectively created idealized microscopic heat spikes.
The researchers ran 20 independent simulations at several energies, including 2, 10, 20, 50, 100 and 200 keV.
Above approximately 30 keV, these artificial heat spikes reproduced the same qualitative transition toward superlinear defect production seen in the full cascade simulations.
That provides an important computational cross-check.
The billion-atom cascade calculations suggest that unusually dense energy deposition is responsible for the superlinear regime.
The artificial experiments isolate that mechanism.
Together, the two approaches strengthen the physical interpretation.
Defects become nanometer-scale structures
The simulations also reveal the physical structures behind the changing damage rate.
The largest defect clusters that frequently form directly in pristine bulk tungsten contain on the order of 1,000 vacancies or self-interstitial atoms.
Depending on morphology, those clusters can span approximately 4–10 nanometers.
At high energies, almost all self-interstitial atoms become part of clusters, while approximately 60% of vacancies are clustered.
Those clusters can include void-like vacancy cores and large dislocation structures.
The computer is therefore doing more than counting defects.
It is revealing the nanoscale architecture of damage.
LUMI becomes part of the physics experiment
The calculations were performed on LUMI, the EuroHPC supercomputer hosted by CSC in Finland.
The GPU implementation of tabGAP was specifically optimized for the accelerator architecture, and the researchers' supplemental data identify the GPU-hours associated with individual cascade simulations.
The work was also carried out partly through the EUROfusion E-TASC Advanced Computing Hub, with access to LUMI awarded by the University of Helsinki.
This is a useful reminder that today's scientific discoveries increasingly depend on the interaction between researchers and computing infrastructure.
The supercomputer isn't simply a place where the calculations happen.
The architecture of the machine shapes which scientific questions can be asked.
The bigger HPC lesson
The most exciting part of this research may ultimately have little to do with tungsten alone.
The work demonstrates a general strategy for computational science:
Use machine learning to make a high-fidelity physical model computationally affordable, optimize it for GPUs, scale the simulation to unprecedented system sizes, and then use the resulting data to discover physics that smaller simulations cannot see.
That strategy is appearing across materials science, fusion, chemistry, astrophysics and climate modeling.
Here, it has pushed molecular dynamics into a regime that was previously inaccessible.
The result isn't simply a faster calculation.
It is a new observation.
From simulation to fusion reactor
Fusion reactor designers ultimately need to know how materials behave over years of neutron exposure.
No computer can simulate every atom in a reactor wall for the reactor's entire operating lifetime.
Instead, scientists need a hierarchy of models.
Atomistic simulations describe the earliest stages of radiation damage.
Those results feed mesoscale models.
Those models inform materials-property predictions.
And those predictions eventually feed reactor-scale simulations and engineering decisions.
The quality of the entire hierarchy depends partly on whether its lowest-level physics is correct.
That is why the new billion-atom calculations matter.
They provide data in a previously inaccessible energy range and reveal that the conventional picture needs to account for four regimes rather than a simple progression from sublinear to linear behavior.
Computing the future of fusion materials
The research characterizes four distinct regimes of radiation-induced damage in tungsten, utilizing billion-atom molecular dynamics simulations to track defect evolution from the displacement threshold to the megaelectronvolt range. By elucidating the transition from sublinear heat-spike effects to superlinear clustering and eventual linear growth at high energies, this study provides a more accurate computational framework for predicting material degradation. These findings enable a refined understanding of tungsten performance, facilitating the optimization of plasma-facing components for the extreme conditions inherent in future fusion reactors.