AutoPas: Auto-Tuning Particle Simulations

• AutoPas is an open-source C++ library for node-level force computation in short-range particle simulations, employing dynamic algorithm tuning.
• It integrates diverse configuration options including containers, traversals, data layouts, and SIMD vectorization to adapt performance across varying simulation conditions.
• Its auto-tuning methods, such as Full Search, Predictive, Expert, and Random Forest, optimize runtime performance based on simulation-specific characteristics.

AutoPas is an open-source C++ library and node-level auto tuning software for efficient short-range particle simulations with a cutoff distance, like MD simulations. Its central design goal is automatic performance portability: it exposes many algorithmic options, including containers, traversals, data layouts, Newton’s third law usage, and, in later work, SIMD vectorization order, and uses online tuning to pick a good combination at runtime. AutoPas is not a full MD package; instead, it focuses on node-level force computation and is typically embedded in higher-level workflows. Although its original design target is short-range molecular dynamics, it has also been used as the core particle back-end of LADDS for deterministic conjunction tracking in long-term space debris simulations (Newcome et al., 6 May 2025, Gall et al., 3 Dec 2025, Gómez et al., 2022).

1. Conceptual model and software role

AutoPas is organized around a black-box particle container abstraction. In this model, simulation code is oblivious to the actual data structure used within AutoPas, which allows AutoPas to choose its internal data structures and algorithms, such as Linked Cells or Verlet Lists, as well as algorithmic optimizations. In short-range MD, the dominant cost per time step is computing pair forces for all particle pairs within a cutoff radius $r_c$; AutoPas addresses this by coupling particle storage, neighbor identification, pairwise iteration, and runtime algorithm selection in one node-level library (Newcome et al., 6 May 2025, Gómez et al., 2022).

The library is designed for shared-memory HPC systems and can be instantiated separately per MPI subdomain. This gives each rank its own AutoPas container and tuning process, so different subdomains can use different algorithmic configurations at the same time. In the MD-focused studies, AutoPas is the force-computation backend rather than the full simulation environment; the higher-level code remains responsible for integrators, MPI domain decomposition, particle migration, ghost exchange, thermostats, and similar workflow logic (Newcome et al., 6 May 2025, Gall et al., 3 Dec 2025).

A recurrent theme across the literature is that there is no single best configuration for all short-range particle simulations. The optimal choice depends on particle density and spatial distribution, cutoff radius and skin, force model, hardware, temporal evolution of the system, and the local properties of each MPI rank. AutoPas was therefore designed around dynamic in-simulation tuning rather than a fixed algorithmic recipe (Newcome et al., 6 May 2025).

2. Containers, traversals, data layouts, and configuration space

AutoPas supports multiple neighbor identification algorithms and traversals. In the algorithm-selection study, the configuration space includes Linked Cells (LC), Verlet Lists (VL), Verlet List Cells, and Pairwise Verlet Lists as particle containers, together with traversals such as List_Iter, C01, C08, C18, C04_HCP, C04, SLI_C02, and SLI. The same study considers Cell Size Factor values $\{0.5, 1\}$ for LC-like containers, AoS and SoA data layouts, and Newton’s third law enabled or disabled, yielding 116 compatible configurations (Newcome et al., 6 May 2025).

The configuration can be written as

$$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$

subject to compatibility constraints. This formulation makes explicit that AutoPas is not a single algorithm but a family of algorithms whose runtime behavior depends on the chosen tuple $\theta$ and on the current simulation state (Newcome et al., 6 May 2025).

The data-layout layer is especially important for SIMD and cache behavior. AoS (Array of Structures) is std::vector<Particle>, with each particle storing position, velocity, force, and related fields. SoA (Structure of Arrays) stores each property in a separate contiguous array, such as $x[i]$, $y[i]$, $z[i]$, and $F_x[i]$, which is better suited for SIMD load/store patterns. In the SIMD-focused work, vectorization assumes SoA because SIMD load/store cost dominates otherwise (Gall et al., 3 Dec 2025).

Later work extends this configuration space one level deeper into the SIMD kernel. AutoPas adds vectorization pattern as another tunable parameter, so the search space becomes

$$[\text{Container}] \times [\text{Traversal}] \times [\text{Newton3}] \times \dots \times [\text{Data Layout}] \times [\text{VectorizationPattern}],$$

with patterns such as $1/v$, $\{0.5, 1\}$0, $\{0.5, 1\}$1, and $\{0.5, 1\}$2 on a $\{0.5, 1\}$3-lane SIMD unit (Gall et al., 3 Dec 2025).

3. Dynamic auto-tuning and algorithm selection

The original AutoPas tuning mechanism is Full Search. At each tuning phase, every compatible configuration is run for a small number of iterations, a performance metric is measured, and the fastest configuration is selected until the next tuning phase. The tuning objective combines container build/update cost and force-calculation cost; forces computed during tuning are not discarded, so the overhead comes from spending iterations in suboptimal configurations rather than from separate benchmarking runs (Newcome et al., 6 May 2025).

Because Full Search may trial many very slow configurations, later work develops three alternative selection strategies: Predictive, Expert, and Random Forest. The Predictive strategy uses recent timing history for each configuration and fits a simple linear model to forecast near-future performance; the Expert strategy uses fuzzy logic and hand-crafted rules derived from live simulation statistics; the Random Forest strategy uses scikit-learn’s RandomForestClassifier trained on synthetic particle distributions and live occupancy statistics. These approaches can achieve speedups of up to 4.05 compared to prior approaches and 1.25 compared to a perfect configuration selection without dynamic algorithm selection (Newcome et al., 6 May 2025).

The evaluation also shows that dynamic selection is not uniformly advantageous. In Heating Sphere, Expert and Random Forest achieve about $\{0.5, 1\}$4 speedup relative to the best static configuration because the simulation transitions from a dense sphere to a sparse gas and the optimal AutoPas configuration changes accordingly. In Exploding Liquid and Rayleigh–Taylor, however, the configuration that dominates the total runtime on the heaviest-loaded ranks remains relatively stable, so dynamic tuning can match or modestly lag the best static configuration because of remaining tuning overhead (Newcome et al., 6 May 2025).

A common misconception is therefore that runtime tuning always improves over the best fixed configuration. The published results do not support that claim. What they support is narrower and more precise: when the optimum changes with time or differs substantially across ranks, dynamic selection can outperform static choices; when a single configuration remains dominant for the most expensive parts of the run, the gains can be small or absent (Newcome et al., 6 May 2025).

4. SIMD vectorization as a tunable dimension

The SIMD work extends AutoPas’s auto-tuning philosophy into the force functor itself. Instead of having a single fixed way of vectorizing the pairwise kernel, AutoPas treats the vectorization order—how many $\{0.5, 1\}$5-particles and $\{0.5, 1\}$6-particles are loaded into SIMD registers—as a tunable parameter that can be switched on the fly based on the current simulation state and the chosen neighbor algorithm. The implementation uses Google’s Highway SIMD wrapper, which provides ISA-independent source code across AVX2, AVX-512, and related instruction sets (Gall et al., 3 Dec 2025).

On an AVX-512 double-precision machine with vector length $\{0.5, 1\}$7, the study evaluates four patterns: 1/VectorLength, VectorLength/1, 2/VectorLengthDiv2, and VectorLengthDiv2/2. These patterns trade off register filling, blank-lane overhead on short lists, reduction cost, scalar versus vector stores, and the handling of Newton3 updates. The results show that the best pattern is highly context-dependent. For Linked Cells, VectorLength/1 can be up to 25% slower than other orders for smaller cutoffs but up to ~12% faster for larger cutoffs; with traversal lc_c08, VectorLength/1 has 22–33% lower runtime than the other patterns; with Newton3 enabled, 1/VectorLength is up to 15% faster than the others (Gall et al., 3 Dec 2025).

For Verlet Cluster Lists, the dominant factor is cluster-size alignment with the vector pattern. With cluster size $\{0.5, 1\}$8 on an 8-way machine, VectorLengthDiv2/2 performs best, up to 45% faster than the others. With cluster size $\{0.5, 1\}$9, patterns aligned to cluster size yield up to 30% speedup over non-aligned patterns. In the multi-rank Exploding Liquid benchmark, different MPI ranks choose different optimal vectorization patterns as local density changes, and tuning for energy rather than runtime can select a different pattern, with energy speedup up to ~2.2× over 1/VectorLength on one rank (Gall et al., 3 Dec 2025).

This work clarifies that AutoPas’s “algorithm configuration” is not limited to containers and traversals. It includes micro-kernel choices inside the force loop, and the optimal SIMD mapping depends on cutoff, traversal, Newton3 usage, cluster size, local density, and the optimization target itself—runtime or energy (Gall et al., 3 Dec 2025).

5. Extension to three-body interactions and multiple time stepping

AutoPas has also been extended beyond pairwise Lennard–Jones kernels. In the three-body study, it serves as the platform for a novel shared-memory three-body cutoff method and for the integration of r-RESPA multiple time stepping. The cutoff-based branch uses AutoPas for both two-body and three-body interactions on a single node with up to 72 hardware threads of a dual-socket Xeon Platinum 8360Y node, whereas a separate DirectSum branch implements a communication-reducing distributed-memory algorithm for no-cutoff three-body interactions (Martin et al., 15 Jul 2025).

For two-body interactions, the relevant potential is the Lennard–Jones $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$0 potential. For three-body interactions, the study implements the Axilrod–Teller–Muto (ATM) potential

$$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$1

with the cutoff condition

$$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$2

AutoPas uses the same linked-cells structure for both cases but adds a three-body C01 traversal that enumerates valid cell-offset patterns $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$3 around a base cell and generates particle triplets on the fly. The paper explicitly states that, in this traversal, Newton’s Third Law of motion cannot be used, because each processor processes a base cell and accesses the surrounding cells in all directions (Martin et al., 15 Jul 2025).

The r-RESPA formulation splits forces into a fast two-body part and a slow three-body part,

$$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$4

with the two-body force evaluated every time step and the three-body force evaluated only every $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$5-th step, where $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$6 is the step-size factor. AutoPas supplies the cutoff-based $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$7 and $$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$8 kernels, while the application-level integration loop schedules them according to the r-RESPA algorithm (Martin et al., 15 Jul 2025).

The measured performance and accuracy trade-offs are specific. In the aluminum scenario, a cutoff-based two-body-plus-three-body simulation requires ~13.4× more runtime than a pure two-body cutoff simulation, whereas the analogous DirectSum case is ~393.6× more expensive than pure two-body direct simulation. Step-size factors from 2 to 6 show little change in RDF and pressure relative to full Störmer–Verlet, while 12 causes clear deterioration in energy behavior and a noticeable shift in pressure. The study therefore positions AutoPas as a research platform for advanced MD algorithms, not only standard two-body short-range interactions (Martin et al., 15 Jul 2025).

6. Cross-domain use, demonstrated performance, and limitations

AutoPas’s use in LADDS shows that it is not restricted to thermal MD. In that architecture, the implementation has two major components: storage for all simulated particles, implemented with AutoPas, and a numerical propagator, implemented in a separate library called OrbitPropagator. AutoPas stores debris and satellite objects, provides neighbor search and pairwise iteration, and is treated as a black-box particle container. The conjunction tracking algorithm is implemented as an AutoPas functor passed to the pairwise interaction interface; AutoPas performs neighbor search within a cutoff distance and applies the functor to candidate pairs, while the higher-level code performs orbit propagation and simulation control (Gómez et al., 2022).

The simulation uses a finite box container because AutoPas is restricted to finite, box-shaped domains. The specific setup is a cube with side length 20000 km centered at Earth, with no boundary wrapping. The conjunction kernel computes the time of closest approach by sub-timestep linear interpolation and tests

$$\theta = (\text{Container}, \text{Traversal}, \text{CSF}, \text{Layout}, \text{N3L}),$$9

so the debris problem is cast as a short-range pairwise interaction analogous to a molecular-dynamics force functor. This suggests a broader methodological point: AutoPas can be used in orbital dynamics and space debris environment modelling as long as interactions are short-range and pairwise (Gómez et al., 2022).

The proof-of-concept simulation propagates 16,024 particles in low-Earth orbit over 20 years with time step $\theta$0 s. The runtime breakdown is explicit: Integrator accounts for >61% of total time, Collision detection for ~20%, Burnup detection for ~15%, Container update for <3%, and I/O and initialization for <1%. Performance is about 0.01 s wall time per simulation iteration, corresponding to ~10.91 hours per simulated year, and the full 20-year run completes in 218.25 hours on 28 cores of a CoolMUC2 node (Gómez et al., 2022).

These results are often summarized as evidence that modern computational tools finally enable deterministic conjunction tracking. In the AutoPas-specific reading, the important point is narrower: the AutoPas pairwise step is efficient enough that, in this application, the integrator rather than the neighbor search becomes the main bottleneck. At the same time, the cited literature is explicit about current limitations: Full Search can incur substantial overhead; expert-rule systems are labor-intensive; Random Forest selection depends on training data and hardware; the three-body C01 traversal does not exploit Newton’s third law; the cited studies do not emphasize GPU support; and the LADDS use case is confined to box-shaped domains because of AutoPas’s current domain restriction (Newcome et al., 6 May 2025, Martin et al., 15 Jul 2025, Gómez et al., 2022).

Across these applications, AutoPas appears less as a single MD code than as a tunable computational substrate for short-range particle problems. Its distinctive contribution is the combination of container abstraction, traversal diversity, data-layout control, shared-memory parallelism, and runtime algorithm selection, extended in later work from macro-level container choice down to the SIMD packing order inside the force kernel (Gall et al., 3 Dec 2025).