TurboGAP: High-Performance GAP & MD Simulations

• TurboGAP is a computational framework that combines high-performance Gaussian Approximation Potential fitting with large-scale molecular dynamics simulation for accurate atomistic modeling.
• It efficiently integrates experimental-agreement forces using a modified molecular dynamics (MAD) approach, enabling simultaneous minimization of potential energy and experimental deviation.
• The framework leverages both MPI+OpenMP and GPU implementations to achieve linear scaling and significant speedups in simulations, fitting, and radiation-damage dynamics.

Searching arXiv for recent TurboGAP papers and related records. TurboGAP denotes a set of closely related high-performance implementations for Gaussian Approximation Potential (GAP) workflows. In the cited literature, the name refers both to a data-parallel implementation of GAP fitting using ScaLAPACK with MPI and OpenMP, and to a software package for efficient molecular dynamics simulations using GAP machine-learning interatomic potentials (MLIP). In its 2025 molecular augmented dynamics (MAD) formulation, TurboGAP evaluates the GAP interatomic energy together with analytical “experimental forces” derived from an experimental-agreement potential, with linear-scaling routines for X-ray/neutron diffraction and local observables such as core-electron binding energies used in X-ray photoelectron spectroscopy. The same code base is also extended to radiation-induced cascade simulations through two-temperature molecular dynamics, adaptive timestepping, and electronic stopping models, with demonstrated simulations from $10^3$ to $10^6$ atoms (Klawohn et al., 2022, Zarrouk et al., 26 Sep 2025, Saha et al., 30 Sep 2025).

1. Scope and problem setting

TurboGAP is positioned at the intersection of MLIP fitting, large-scale molecular dynamics, and experiment-driven atomistic modeling. For disordered systems, the central difficulty is that aligning theoretical atomistic structural models with available experimental data requires navigating a vast configurational space, while faithful realizations require large system sizes with quantum-mechanical accuracy in order to capture the distribution of structural motifs present in experiment. Traditional equilibrium sampling approaches offer no guarantee of generating structures that coincide with experimental data for such systems. Within that context, MAD is introduced as a modified molecular dynamics method that can generate ab-initio accurate, low-energy structures through a multi-objective optimization of the interatomic potential energy and the experimental potential, with TurboGAP acting as the high-performance engine for the required force evaluations (Zarrouk et al., 26 Sep 2025).

A distinct, but complementary, use of the same name appears in the GAP-fitting literature. There, TurboGAP is a data-parallel implementation of Gaussian Approximation Potential fitting that overcomes the single-node memory and compute limits of earlier serial or purely OpenMP-parallel codes, scales to thousands of cores, and enables higher-level workflows such as committee models or hyperparameter optimisation (Klawohn et al., 2022).

A common misunderstanding is to treat TurboGAP as only a force-evaluation backend or only a training code. The published record instead assigns it both roles: distributed fitting of GAP models, and large-scale deployment of those models in molecular dynamics, including MAD and radiation-damage simulations. This suggests that the term “TurboGAP” is best understood as a computational framework centered on GAP models rather than as a single narrowly scoped executable.

2. GAP formalism and interatomic representation

At the level of the interatomic model, TurboGAP uses the standard GAP decomposition in which the total potential energy is written as a sum of atomic energies,

$E_{\rm tot}=\sum_i E_i,$

or, equivalently in the fitting formulation,

$E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$

Each local contribution is obtained through sparse Gaussian-process regression over a set of reference environments. In the simulation-oriented description,

$$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$

where $\mathbf{D}_i$ is the vector of descriptors, $K$ is the kernel, $\{\alpha_s\}$ are regression coefficients, and $N_{\rm sparse}$ is the number of sparse points. In the fitting-oriented description, the sparse model is written as

$\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$

with a sparse set $10^6$0 and weights $10^6$1 determined from energies, forces, and stresses through a regularised least-squares formulation (Saha et al., 30 Sep 2025, Klawohn et al., 2022).

The silicon GAP used for radiation-damage studies combines three separate descriptors: a two-body kernel up to $10^6$2 Å, a three-body kernel up to $10^6$3 Å, and a turboSOAP many-body kernel with hard cut $10^6$4 Å, smooth cutoff $10^6$5 Å, $10^6$6, $10^6$7, and $10^6$8. To capture the highly repulsive part of the potential-energy surface at very short distances, a tabulated Ziegler–Biersack–Littmark (ZBL) potential is appended,

$10^6$9

with the universal screening function

$E_{\rm tot}=\sum_i E_i,$0

This combination is specifically motivated by the need to describe cascade events involving very short interatomic separations (Saha et al., 30 Sep 2025).

In the fitting formulation, numerical stability is improved by constructing an over-determined Tikhonov-regularised least-squares system and solving it through distributed QR factorisation rather than by direct inversion of the normal equations. That choice is central to TurboGAP’s ability to handle large sparse sets and large observation counts (Klawohn et al., 2022).

3. Molecular augmented dynamics and experimental-force machinery

In MAD, TurboGAP couples the GAP potential to an experimental-agreement potential within a single Hamiltonian,

$E_{\rm tot}=\sum_i E_i,$1

where $E_{\rm tot}=\sum_i E_i,$2 is the GAP interatomic energy and

$E_{\rm tot}=\sum_i E_i,$3

is a quadratic experimental potential measuring deviation of $E_{\rm tot}=\sum_i E_i,$4 predicted observables from experiment. The MAD dynamics then minimizes both the physical energy and the experimental deviation simultaneously. TurboGAP provides both the GAP forces $E_{\rm tot}=\sum_i E_i,$5 and the analytical gradients, or “experimental forces,” $E_{\rm tot}=\sum_i E_i,$6 in $E_{\rm tot}=\sum_i E_i,$7 time, enabling large-scale simulations up to millions of atoms (Zarrouk et al., 26 Sep 2025).

The linear-scaling character of the method depends on the fact that all key routines in TurboGAP—GAP evaluation, neighbor-list construction, and experimental observable gradients—scale linearly with system size $E_{\rm tot}=\sum_i E_i,$8. Neighbor searching with a finite cutoff $E_{\rm tot}=\sum_i E_i,$9 is $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$0. Kernel-density estimation of partial pair distributions and their gradients for XRD, ND, and PDF is $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$1. Summations over local descriptors for XPS are $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$2. The overall cost is therefore described as

$E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$3

The CPU algorithmic steps per MD step are specified explicitly: build or update the neighbor list; compute $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$4 via Gaussian KDE; Fourier transform partial PDFs to structure factors $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$5 via a short-cut explicit sum over $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$6-points; assemble predicted intensity $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$7 and its gradient; combine with experimental weights and form $E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$8; then return total forces and virial contributions. This sequence makes clear that experimental matching is not a post-processing layer but part of the force loop itself (Zarrouk et al., 26 Sep 2025).

TurboGAP also implements the generalized virial stress including experimental forces,

$E_A=\sum_{i=1}^{N_A}\varepsilon(x_i).$9

$$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$0

Because $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$1 gathers contributions from both $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$2 and $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$3, TurboGAP can run fully flexible-cell NPT MAD simulations in which volume and shape evolve to match experimental densities automatically (Zarrouk et al., 26 Sep 2025).

4. CPU and GPU implementations and measured performance

On CPUs, TurboGAP uses an MPI + OpenMP programming model. Each MPI rank owns a subset of atoms; within each rank, OpenMP threads loop over local atoms. Neighbor lists are Verlet lists with cutoff $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$4, updated infrequently relative to MD timesteps. Data structures include SoA arrays of positions $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$5, species labels, partial-PDF histograms $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$6 in radial bins, and pre-computed scattering factors $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$7 or neutron lengths $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$8. Strong-scaling on 27 000 atoms achieves $$E_i=\sum_{s=1}^{N_{\rm sparse}} \alpha_s\,K(\mathbf{D}_i,\mathbf{D}_s),$$9 efficiency up to 512 CPU cores, and the experimental-force subroutines remain strictly $\mathbf{D}_i$0 and comparably priced even when the dominant cost remains the GAP force and descriptor evaluation (Zarrouk et al., 26 Sep 2025).

The GPU path is implemented in HIP, with HOP for runtime portability, plus cuBLAS/hipBLAS for small dense linear algebra. Positions and neighbor lists are transferred to GPU once, while partial-PDF accumulators reside in GPU global memory. For very large $\mathbf{D}_i$1, the domain of atoms is split into $\mathbf{D}_i$2 batches to fit HBM memory, with typical $\mathbf{D}_i$3. The GPU workflow is decomposed into three kernels: PDF build; PDF-to-structure-factor transform through a small dense matrix–vector product; and intensity-plus-gradient contraction. Host–device transfer is overlapped with compute (Zarrouk et al., 26 Sep 2025).

Benchmarks on the LUMI supercomputer quantify both absolute timing and scaling. For strong scaling on 27 000 atoms, the XRD/ND experimental-force subroutines decrease from $\mathbf{D}_i$4 s on 1 rank ($\mathbf{D}_i$5 cores) to $\mathbf{D}_i$6 s on 512 ranks on the CPU side, and from $\mathbf{D}_i$7 s on 1 GPU to $\mathbf{D}_i$8 s on 4 GPUs on the GPU side. Relative to 1 $\mathbf{D}_i$9 64-core CPU, the reported speedups are $K$0 for 512 CPU cores and $K$1 for 4 GPUs. For weak scaling from 1 728 to $K$2 atoms, both CPU and GPU timings remain linear, with $K$3 CPU efficiency and $K$4 GPU efficiency at 512 ranks. On 1 $K$5 AMD MI250x die, the experimental-force subroutines for $K$6 run in a few tens of milliseconds—up to $K$7 faster than a single 64-core MPI rank on LUMI-C. In a full MAD MD step, adding the GPU experimental forces costs $K$8 of the total, versus $K$9 on CPU, and the GPU version runs $\{\alpha_s\}$0 faster per MD step than the CPU one, including GAP evaluation (Zarrouk et al., 26 Sep 2025).

A frequent misconception is that the experimental observables dominate the runtime once they are coupled to a machine-learned potential. The published measurements indicate the opposite for the GPU implementation: explicit batching and high-throughput kernels push the experimental-force cost to $\{\alpha_s\}$1 of the total GAP+MD expense (Zarrouk et al., 26 Sep 2025).

5. Radiation-induced cascade simulations

For radiation damage simulations, TurboGAP is extended with a two-temperature molecular dynamics model, also denoted TT-MD or “EPH,” based on electron density-dependent coupling of electronic and atomic subsystems. The atomic subsystem follows a generalized Langevin equation,

$\{\alpha_s\}$2

while the electronic subsystem obeys

$\{\alpha_s\}$3

In TurboGAP these equations are integrated jointly via finite-difference for the electronic mesh, Verlet, and stochastic updates (Saha et al., 30 Sep 2025).

Two electronic stopping models are implemented. The friction-based electronic stopping (FES) model applies a velocity-dependent drag using a tabulated stopping power $\{\alpha_s\}$4, with a user-selectable low-energy cutoff $\{\alpha_s\}$5. The first-principles-derived stopping model instead uses the electron-density-dependent friction tensor $\{\alpha_s\}$6, derived from time-dependent DFT studies of stopping power; no ad hoc cutoff is needed, and the same model also produces the correct electron–phonon coupling at low energies. TurboGAP additionally implements adaptive calculation of the timestep through maximum displacement and maximum energy-change criteria, with the actual $\{\alpha_s\}$7 clamped to $\{\alpha_s\}$8 (Saha et al., 30 Sep 2025).

The published silicon cascade setup uses PKA energies of 0.1, 0.2, 0.4, 1, 2, 5, and 10 keV, with system sizes from 4 096 atoms up to $\{\alpha_s\}$9 atoms. Initial equilibration is at 300 K and 0 bar for 20 ps with a Berendsen thermostat/barostat. The EPH model employs an electronic bath approximately $N_{\rm sparse}$0 larger than the atomistic cell with 20–22 Å grid spacing in $N_{\rm sparse}$1, $N_{\rm sparse}$2, and $N_{\rm sparse}$3. FES runs use an SRIM-2013 $N_{\rm sparse}$4 file and two cutoffs, $N_{\rm sparse}$5 eV. Production cascades run to 10 ps (Saha et al., 30 Sep 2025).

Defect generation is analyzed by Wigner–Seitz analysis in OVITO. The number of interstitials peaks during the “heat spike” at 0.2–0.7 ps and then relaxes to survivors. EPH yields systematically more survivors, for example $N_{\rm sparse}$6 versus $N_{\rm sparse}$7 for 2 keV, and larger average clusters. Interstitial clusters are defined by a 5.3 Å criterion; violin plots show that EPH’s largest clusters exceed 100 atoms, whereas FES stays below $N_{\rm sparse}$8. Ion-beam mixing is quantified through

$N_{\rm sparse}$9

with a spectrum-integrated version for heavy-ion experiments. The reported $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$0 values agree within error bars with experiments of Matteson et al. (1981) and Paine et al. (1981); for Ar 110 keV, $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$1 $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$2 versus $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$3 $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$4 (Saha et al., 30 Sep 2025).

6. Massively parallel fitting and workflow integration

The fitting implementation of TurboGAP distributes the training structures and their associated observations almost evenly across MPI ranks and uses OpenMP within each rank to parallelise over sparse points. Data partitioning is along the observation dimension $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$5, since both memory and computation scale proportionally to $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$6. A greedy load-balancing assigns each XYZ structure to the MPI rank currently holding the fewest total observations. Each rank then computes the local descriptor vectors and their kernels for all sparse points, forms its partial covariance block, and assembles the local contribution to the regularised least-squares system without any inter-rank communication during descriptor or kernel evaluation (Klawohn et al., 2022).

The global solve is mapped onto ScaLAPACK. TurboGAP distributes the rows of the least-squares matrix $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$7 and the entries of $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$8 in a block-cyclic fashion over a $\varepsilon(x_i)=\sum_{j=1}^{M} c_j\,k(x_i,x_j),$9 MPI grid; applies PDGEQRF for QR factorisation; applies PDORMQR to form $10^6$00; and solves the triangular system through PDTRTRS. To minimise memory blow-up in PDTRTRS, the code dynamically switches the ScaLAPACK column block size for that call. A practical column block size of order 50–200 is reported to provide a good trade-off between workspace and performance (Klawohn et al., 2022).

The benchmarked training sets illustrate the scale of the distributed fit. For a High-Entropy Alloy MoNbTaVW dataset with $10^6$01 observations and $10^6$02 sparse points, the global matrix $10^6$03 has approximately $10^6$04 entries, or about 67.6 GB. For a Silicon Carbide dataset with $10^6$05 observations and the same $10^6$06, $10^6$07 has approximately $10^6$08 entries, or about 406 GB. The covariance time $10^6$09 scales almost perfectly inversely with the total number of cores, whereas the QR solve time $10^6$10 benefits less from node count because of inter-node communication overheads in ScaLAPACK. A 50:50 split of MPI tasks to OpenMP threads per node generally maximises efficiency and keeps memory overhead low, while pure MPI often saturates at moderate node counts and incurs a large memory footprint (Klawohn et al., 2022).

The reduction in fit times and the removal of single-node memory limits are explicitly connected to higher-level workflows. The cited applications are committee-model uncertainty quantification, automated hyperparameter optimisation, and on-the-fly active learning loops combining DFT sampling with GAP retraining. Best-practice recommendations include performing sparse-point selection in a small serial pre-run, using a moderate ScaLAPACK column block size of around 100, and reserving some cores for OpenMP to amortise memory and communication cost (Klawohn et al., 2022).