Improved capabilities of the TurboGAP code for radiation induced cascade simulations: an illustration with silicon

Technical Summary: Improved Capabilities of the TurboGAP Code for Radiation Induced Cascade Simulations

Problem Statement
Accurate atomistic simulation of radiation damage requires balancing computational efficiency with the physical fidelity of the potential energy surface (PES). While Machine Learning Interatomic Potentials (MLIPs), such as the Gaussian Approximation Potential (GAP), offer superior accuracy compared to empirical potentials, their application to radiation cascade simulations is often limited by computational cost, particularly for high primary knock-on atom (PKA) energies. Furthermore, standard Molecular Dynamics (MD) codes often lack essential features required for realistic radiation damage modeling: the ability to handle adaptive timesteps during high-energy collisions, effective cooling of simulation cell boundaries, and accurate modeling of electronic energy dissipation (electronic stopping). Traditional friction-based electronic stopping (FES) models often rely on arbitrary kinetic energy cutoffs, which can significantly influence defect microstructures. There is a need for a software framework that combines the efficiency of TurboGAP with advanced dissipation models, such as the electron-phonon (EPH) coupling model, to simulate cascades in systems containing up to $10^6$ atoms with PKA energies up to 10 keV.

Methodology
The authors enhanced the TurboGAP software package by implementing three critical modules:

1. Electronic Energy Dissipation Models:

   • Friction-based Electronic Stopping (FES): Implemented as a velocity-dependent friction force derived from pre-calculated electronic stopping power (ESP) tables (e.g., from SRIM). A user-defined kinetic energy cutoff ($E_{el\_cut}$) determines which atoms experience this force.
   • Electron-Phonon (EPH) Model: Integrated into a generalized two-temperature molecular dynamics (TTMD) framework. This model couples the atomic and electronic subsystems via electron density-dependent parameters. It solves the heat diffusion equation for the electronic temperature ($T_e$) simultaneously with atom dynamics. The atomic motion is governed by adiabatic forces, friction forces dependent on relative velocities, and stochastic random forces derived from the fluctuation-dissipation theorem. This approach avoids arbitrary kinetic energy cutoffs.

2. Adaptive Timestep Calculation: To capture the steep potential energy gradients during close collisions, the timestep is dynamically recalculated based on two criteria: the maximum allowed displacement of any atom ($x_{max}$) and the maximum change in kinetic energy ($e_{max}$). The minimum of the timesteps derived from these criteria is selected.

3. Atom Grouping and Border Cooling: A grouping mechanism allows for the selection of atoms based on spatial criteria (e.g., a border layer) to apply thermostats. This enables the dissipation of excess heat from the cascade into the bulk material, maintaining the desired simulation temperature.

Potential Development
A specific GAP MLIP for Silicon (Si) was retrained using 2-body, 3-body, and turboSOAP descriptors (an optimized variant of the SOAP descriptor for improved efficiency). The potential was trained on a dataset including various Si structures, with specific regularization parameters for liquid and amorphous phases. To accurately model high-energy repulsive interactions during collisions, the universal Ziegler-Biersack-Littmarck (ZBL) potential was appended to the GAP model as an external potential, active at short interatomic distances.

Key Results

• Performance and Scaling: The new modules were benchmarked on High-Performance Computing (HPC) setups. Scaling tests with up to 600 processors for systems of $10^6$ atoms demonstrated that the implemented modules maintain the computational efficiency of the original TurboGAP code.
• Code Validation: Comparisons with LAMMPS (using the USER-EPH plugin and standard GAP implementations) showed excellent agreement in the evolution of thermodynamic properties (temperature, pressure, energy) and cumulative electronic energy loss for 100 eV PKA cascades.
• Defect Production and Clustering (2 keV PKA):
  • The EPH model predicted a higher number of surviving defects (33) compared to FES models with 10 eV (26) and 1 eV (~22) cutoffs.
  • The EPH model resulted in larger defect clusters (some exceeding 100 interstitials), whereas FES models produced clusters generally below 100 interstitials.
  • The FES model's results were sensitive to the arbitrary choice of kinetic energy cutoff, while the EPH model results were robust across different electronic subsystem settings.
  • The EPH model allowed for continuous energy transfer to the electronic system even after the cascade cooled, unlike the FES model where energy transfer ceased once atoms dropped below the cutoff energy.
• Ion-Beam Mixing:
  • Mixing values ($Q$) calculated using the EPH model and GAP MLIP showed good agreement with experimental data for Si irradiated with Ne, Ar, Kr, and Xe ions.
  • The GAP MLIP predictions differed significantly (by a factor of 2–4) from those using the classical Stillinger-Weber (SW) potential, highlighting the importance of accurate PES representation.
  • The EPH model's mixing values were comparable to FES results but achieved without the arbitrariness of kinetic energy cutoffs.

Significance and Claims
The paper claims that the integration of the EPH model, adaptive timesteps, and border cooling into TurboGAP significantly enhances the capability for radiation damage simulations. By combining the computational efficiency of TurboGAP with the accuracy of GAP MLIPs, the authors successfully simulated cascades with PKA energies up to 10 keV in systems of up to 1 million atoms.

The authors emphasize that the EPH model offers a distinct advantage over the traditional FES model by eliminating the need for arbitrary kinetic energy cutoffs, which are known to affect defect microstructures. The study demonstrates that realistic electronic stopping models, when coupled with accurate MLIPs, improve the prediction of experimental ion-beam mixing values and provide a more physically consistent description of defect clustering and survival in silicon. The work establishes TurboGAP as a robust tool for multi-scale modeling of radiation-induced damage in materials.