D2.5 - Technical report on performance-portable electrostatics

Technical report on performance-portable electrostatics MultiXscale Deliverable 2.5 Deliverable Type: Report Delivered in June, 2025 MultiXscale EuroHPC Centre of Excellence for Multiscale Modelling Acknowledgement Funded by the European Union. This work has received funding from the European High Performance Computing Joint Undertaking (JU) under grant agreement No 101093169. Disclaimer Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European High Performance Computing Joint Undertaking (JU). Neither the European Union nor the granting authority can be held responsible for them. MultiXscale Deliverable 2.5 Page ii Project and Deliverable Information Project Title MultiXscale: EuroHPC Centre of Excellence for Multiscale Modelling Project Ref. Grant Agreement 101093169 Project Website https://www.multixscale.eu EuroHPC Project Officer Matteo Mascagni Deliverable ID D2.5 Deliverable Nature Report Dissemination Level Public Contractual Date of Delivery Project Month 30 (30th June, 2025) Actual Date of Delivery 27th June, 2025 Description of Deliverable Review of the development of performance-portable electrostatic solvers developed as a community library. Document Control Information Document Title: Technical report on performance-portable electrostatics ID: D2.5 Version: As of June, 2025 Status: Accepted by Steering Committee Available at: https://www.multixscale.eu/deliverables Document history: Internal Project Management Link Review Review Status: Reviewed Authorship Written by: Rodrigo Bartolomeu (FZJ) Contributors: Godehard Sutmann (FZJ) Reviewed by: Alan Ó Cais (UB) Approved by: Alan Ó Cais (UB) Document Keywords Keywords: MultiXscale, HPC, Performance Portability, Electrostatic, Load Balance, ALL 27th June, 2025 Disclaimer: This deliverable has been prepared by the responsible Work Package of the Project in accordance with the Consortium Agreement and the Grant Agreement. It solely reflects the opinion of the parties to such agreements on a collective basis in the context of the Project and to the extent foreseen in such agreements. Copyright notices: This deliverable was co-ordinated by Rodrigo Bartolomeu1(FZJ) on behalf of the MultiXscale consortium with contributions from Godehard Sutmann (FZJ). This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit: http://creativecommons.org/licenses/by/4.0 cb 1r[email protected] MultiXscale Deliverable 2.5 Page iii Contents Executive Summary 1 1 Introduction 2 1.1 Scope of the deliverable ................................................ 2 1.2 Target audience . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Report outline ...................................................... 2 1.4 Partner contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Performance Portability 3 3 Electrostatics Methods and Implementation Considerations 7 3.1 Ewald Summation Method .............................................. 7 3.2 Particle-Particle Particle-Mesh Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.3 Implementation Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4 Electrostatic Solver Library 10 4.1 Library Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.2 MPI Communication management . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.3 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5 Integration of the ALL load balancing library into LAMMPS 14 5.1 ALL interface with LAMMPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.2 General LAMMPS Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.3 Benchmark with balancing methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 6 Conclusion and Perspectives 18 Acknowledgements 19 References 19 List of Figures 1 Sequence of operations and program flow within the N-body problem. The main phase of the program is repeated for Ntiterations, each computing a discrete time step of size dt. P.I.: particle initialization; F.C.: force computation; ∆r : propagation of positions in Velocity Verlet; ∆v : propagation of half step velocities in Velocity Verlet integrator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Interactions per second vs number of particles Non the investigated architectures by the different implementations. Note the different scale of the y-axis for GH200 (b). We do not include the results for OpenACC on the MI 250 since they are about a factor of 100 lower than those from the other programming models and nearly independent of system size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Subdivision of the global MPI communicator into sub-communicators to be used for the computation of the real(r)- and k-space contributions, showing the necessary data transfer between global communicator and each subcommunicator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4 Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 16 million randomly placed particles in a periodic system of size 1003on Juwels Booster (Forschungszentrum Jülich). The resources are evenly split between the r-space and k-space computation. . . . . . . . . 12 5 Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 16 million randomly placed particles in a periodic system of size 1003on Juwels Booster (Forschungszentrum Jülich). Three of four GPUs are employed for the computation of the k-space contribution, the other GPUs are used for the r-space contribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 6 Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 32 million randomly placed particles in a periodic system of size 1263on Juwels Booster (Forschungszentrum Jülich). The resources are evenly split between the r-space and k-space computation. . . . . . . . . 13 7 Single node weak scaling of LAMMPS using JEDI and JUWELS Booster (Forschungzentrum Jülich) for the LJ and Rhodopsin benchmarks in LAMMPS using Kokkos. Scaled LJ simulation containing 55M atoms, and Scaled Rhodopsin simulation containing 16M atoms both with 1 MPI process per GPU. No results are presented for Juwels Booster with 1 and 2 GPUs for the Rhodopsin benchmark due to memory limitation to accommodate the data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 8 Strong scaling of LAMMPS using JEDI for the LJ and Rhodopsin benchmarks in LAMMPS using Kokkos. 16 MultiXscale Deliverable 2.5 Page iv 9 Combined strong and weak scaling of LAMMPS at the pre-exascale supercomputer JEDI (Forschungszentrum Jülich). Solid bars use the straggered global load balancing method from the ALL library integrated into LAMMPS and hatched bars use the recursive bisection load balancing method native to LAMMPS. . 17 List of Tables 1 Overview of implementations on investigated architectures. ✓indicates that benchmarks were performed. ✗indicates that the implementation was not supported on that architecture at the time of the benchmark. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Fastest measured performance data across different system sizes, measured in Giga-interactions per second (GInt s) on each platform for each programming model. The fastest overall measured performance is used as reference value for the application efficiency (%) and is indicated in bold. % (cmp. Eq. 5) is the weighted mean of % over all architectures, assuming % =0, when benchmarks could not be performed. The underlined result is the best achieved value for %. (s) indicates variants using shared memory, while (nd) indicates the SYCL variant using nd-ranges. The last column shows the performance performance portability metric P ¯ P ¯.................................................. 6 MultiXscale Deliverable 2.5 Page 1 Executive Summary With the fast changing High Performance Computing (HPC) landscape, the need of performance portable implementations has grown significantly in the past years. Among the TOP 50 machines on the TOP500 list, only 10% do not rely on Graphical Processing Unit (GPU)s as accelerators for its performance. Moreover, accelerators from multiple vendors hold the top spots on the list. The same picture is also present on European High Performance Computing Joint Undertaking (EuroHPC) and Gauss Supercomputing Centre (GCS) machines. Both have supercomputers with AMD and NVIDIA GPUs of different generations, architectures, and compute capabilities. This striking reality poses a growing demand for performance portability of scientific codes, libraries, and applications. For many domain applications, the efficient use of CPUs alone is not sufficient to extract the available performance and energy efficiency available on modern architectures. As an example, the worlds most energy efficient supercomputer JEDI (EuroHPC/FZJ) is based on Nvidia’s Grace-Hopper Superchip, which combines CPU and GPU. In this system the GPUs can be programmed using Compute Unified Device Architecture (CUDA), however doing so the code would not be able to run on other vendor’s platforms. Achieving performance portability can be accomplished by pairing devices and leveraging offloading of compute kernels through frameworks such as OpenMP, OpenACC, Kokkos, or Raja. Our preliminary research on a performance portable Ewald summation method enlightens the necessity for further comparison between different software frameworks to determine the most effective approach. Therefore, a performance portability study was conducted using the N-body problem as a test case to evaluate the efficacy of various programming models across major GPU vendors. This study aimed to evaluate the relative strengths and weaknesses of each model, providing valuable insights into the optimal approach for achieving performance portability across diverse hardware architectures. For electrostatic computations based on variants of the Ewald summation method, which are often one of the most compute intense and time consuming steps of Molecular Dynamics (MD) simulations, not only portability is desired but also modularity. Task division can be achieved by separating the short-range and long-range contributions to the electrostatic potential and forces. Since the real-space and Fourier k-space parts have different performance characteristics, a multi-domain partitioning approach can be applied in conjunction with adjusting relevant Ewald method parameters. This allows for a balanced computational load across independently running partitions. To achieve optimal load balancing between the real-space and k-space partitions, it is essential to optimize operations performed on both parts as well as the amount of compute resources allocated to each task. In this deliverable, we present advances, implementation details, rationale, and benchmarks related to the aforementioned methods. The outcomes of the performance portability study, as well as the previously reported FFT benchmark, serve as foundation for the design choices included in the development of P3M method and the portable electrostatic solver library. MultiXscale Deliverable 2.5 Page 2 1 Introduction 1.1 Scope of the deliverable This deliverable aims at providing information on the implementation of a performance portable electrostatics solver and the integration of the ALL load balancing library into LAMMPS. This report contains elements relevant to Task 2.3 - Performance-portable electrostatics solvers and load-balancing lead by Forschungszentrum Jülich. 1.2 Target audience This report is intended for readers with technical expertise in HPC or in the use of multi scale particle based simulation software. While potential users might be interested into the details reported here, the main purpose is to report on the results and implementations that might be more of interest for scientific software developers willing to interface their code to the portable electrostatic solver library. 1.3 Report outline This deliverable reports on all relevant topics with respect to Task 2.3, led by Forschungzentrum Jülich. The deliverable is divided into four parts: •Performance Portability: in this section we present the up-to-date status on the performance portability of the most commonly used C++ portable programming models. During the course of Task 2.3 a study was conducted to access the usability and performance of such models. As outcome, an updated status on the portability was reported. Practical experience in the usage as well as a deeper understanding of the implementation details of each model was acquired. This section contributed to the choice of the programming model used for the implementations of electrostatic solvers. •Theoretical foundations of Ewald sum and P3M methods: in this section we present the formulation of the basic Ewald sum, which is the original method to compute the electrostatic energy of an atomistic system under periodic boundary conditions. The accelerated version, P3M, makes use of the Fast Fourier transform, evaluated on a regular mesh, carrying the mapped particle charges. Cardinal B-Splines are applied to map charges to the mesh and to retrieve potential energies and forces back to the particles. While being able to control the approximation error, the P3M method is reducing the computational complexity from O(N3/2) down to O(Nlog N), where Nis the number of particles in the system. •Electrostatic Solver Library: in this section we discuss the design and implementation of the performance portable and scalable library approach for electrostatic solvers. In the present project we focus on the Ewald summation as a benchmark implementation and the P3M as an optimized implementation for production runs on modern HPC systems. Performance portability has been considered for both the real and reciprocal space part of the method. Proper parameter choice for error control and aspects of different scalability of real and reciprocal space parts are discussed and considered in the design of the library. First benchmarks on different architectures are shown and properties of scalability are discussed in terms of different tasks in the method. •Load Balancing: with the intent to benchmark the load balancing methods available in the A Load balancing Library (ALL) within using our key application codes. We have integrated ALL to LAMMPS through an interface that provide users the capability of using extended fix commands. In the context of this deliverable we present benchmarks comparing LAMMPS performance using A100 and GH200 GPUs. 1.4 Partner contributions JSC contributed as planned to the work in this deliverable which is intended to be contribute to realisation of other tasks in the project. MultiXscale Deliverable 2.5 Page 3 2 Performance Portability The ongoing development of new GPUs from various vendors requires continuous monitoring of portable frameworks’ performance and versatility. A primary concern is whether a portable framework is available and functional on a given architecture. While the concept of portability between different platforms is straightforward, defining and measuring performance portability is more complex. Successful program execution is a clear indicator of portability, but evaluating performance portability requires a more nuanced approach. A common method for assessing the performance and capabilities of portable programming models involves using mini-apps or full applications to gather relevant data [1, 2]. This can be achieved by porting a single application to various performance-portable programming models [3–6] or by porting multiple mini-apps to one performanceportable programming model and comparing the results to other available versions [7]. In our study, we focus on the N-body problem [8]. The N-body problem is a significant problem class and has been identified as one of the seven original dwarfs in Ref. [9], each representing a compelling use case with its own challenges. Many scientific fields require computing forces between pairs of interacting particles, such as gravitational forces in astrophysics, electrostatic forces in charged systems, and van der Waals forces for modeling atomic systems. The gravitational potential between particles is given by: U(di j )=−Gmimj di j (1) where G=6.6743·10−11 m3 kg·s2is the gravitational constant, di j =|ri j |is the distance between particles iand j,ri j = ri−rj, and miand mjare the masses of the respective particles. Given a potential U(di j ), the force fi j between particles iand jis given by Fi j =−∇U(di j )=−Gmimj d3 i j ri j . (2) The total force acting on particle iis then given by the sum over all other particles Fi=∑︂ j=i Fi j . (3) To obtain the trajectory of each particle, Newton’s equations of motion x ˙i=vi,v ˙i=ai=Fi mi , (4) are integrated, where the dot-notation represents the time derivative, aiis the acceleration, and miis the mass of particle i. The gravitational N-body problem (the 4th dwarf in Ref. [9]) is a representative of an algorithmic method that includes long-range interactions in open-boundaries and involves summation over all particle pairs in the system to compute the forces on each particle, resulting in a natural computational complexity of O(N2). To avoid this quadratic complexity, more efficient methods like the Barnes-Hut tree method or the fast multipole method (FMM) can be used, reducing the complexity to O(Nlog(N)) or O(N). However, since these methods have a high level of complexity and their performance strongly depends on the implementation and optimization, we focus here on a simple implementation based on the all-pairs computation of energies and forces, which isolates the performance outcome from algorithmic issues. To evaluate the utilization of the hardware architectures in our benchmarks, we compared the throughput of particleparticle computations, where the number of interactions scales as N2. The study used the N-body problem as a test from compute intense algorithm that is present in many MD software. An illustration of the program workflow is presented in Figure 1. In order to evaluate the performance portability of different programming models, the discussed algorithm was implemented in eleven variations. The C++ frameworks considered were: Kokkos [10], OpenMP, OpenACC, SYCL, HIP, CUDA, and pSTL. The implementations were tested on four platforms. An overview of the tested platforms is presented in Table 1. For details about platform specifications, compilers, implementations, and optimizations, readers may refer to [8]. MultiXscale Deliverable 2.5 Page 4 init main I/O Δ𝑣, Δ𝑟 Start P.I. F.C. End Init. F.C. Δ𝑣 𝑁𝑡time steps Figure 1: Sequence of operations and program flow within the N-body problem. The main phase of the program is repeated for Ntiterations, each computing a discrete time step of size dt. P.I.: particle initialization; F.C.: force computation; ∆r : propagation of positions in Velocity Verlet; ∆v : propagation of half step velocities in Velocity Verlet integrator. Architecture OpenMP pSTL Kokkos CUDA HIP SYCL OpenACC Nvidia A100 ✓ ✓ ✓ ✓ ✓ ✓ ✓ Nvidia GH200 ✓ ✓ ✓ ✓ ✓ ✓ ✓ AMD MI250 ✓ ✓ ✓ ✗ ✓ ✓ ✓ Intel GPU Max 1100 ✓ ✓ ✓ ✗ ✗ ✓ ✗ Table 1: Overview of implementations on investigated architectures. ✓indicates that benchmarks were performed. ✗indicates that the implementation was not supported on that architecture at the time of the benchmark. As shown in Figure 2, more powerful GPUs require larger system sizes to saturate the processor, indicated by the flattening of the lines in the plots. The results show that variants using shared memory generally perform better than those without it, particularly on the GH200. However, there are some exceptions that require further investigation. Notably, CUDA on the GH200 with shared memory achieved the best results, while SYCL outperformed HIP on the MI250. In contrast, OpenMP and pSTL outperformed SYCL on the Intel platform. Furthermore, OpenACC, OpenMP, and pSTL surprisingly achieved higher performance on the A100 card compared to CUDA with shared memory. These findings align with a study by [11], which found that SYCL’s software stack improvements led to outperforming native programming models on compute-bound kernels. These results are also presented numerically in Figure 2for discussion: Using the mean application efficiency % (see Equation 5), it can be seen that the models designed to be portable perform better than those targeted at specific platforms (HIP, CUDA). Although this is to be expected, it is noteworthy that pSTL and SYCL (using shared memory) show the best portability quality. Kokkos shows comparatively lower performance on the AMD and Intel platforms. Regarding AMD, our result differs from recent studies, where Kokkos showed excellent performance on AMD-based systems [2,12]. This may be due to the variant chosen for its performance on the GH200. %= %A100 +%GH200 2+%MI250 +%Max1100 3. (5) The respective results for P ¯ P ¯(Table 2) show little deviation from the results for % for the portable programming model, but make it clear that CUDA, HIP, and OpenACC are not portable since they are not applicable (NA) per definition for programming models that cannot be executed on one or more of the investigated architectures. The deviation in the percentages for the portable programming models comes from the different ways of averaging (see. Equation 5). For a more detailed analysis of performance portability the number of taken measurements as well as the number of targeted platforms would need to be increased. Despite this, the performed analysis on the presented data already indicates that using an optimized version of a code on a given architecture (in this case the GH200) might not be optimally suited for other architectures, but still performs satisfactory well on the most modern system from all major HPC GPU vendors. Our results show that OpenMP, pSTL, and SYCL are well suited for portable implementations of the N-body problem across various GPU vendors, while Kokkos showed slightly reduced performance on the non-Nvidia platforms. Kokkos MultiXscale Deliverable 2.5 Page 5 CUDA CUDA (s) HIP HIP (s) Kokkos pSTL OpenMP OpenACC SYCL SYCL (s) SYCL-range 3236439631283 50 100 150 200 𝑁 GInt / s (a) A100 3236439631283 100 200 300 400 𝑁 GInt / s (b) GH200 3236439631283 50 100 150 200 𝑁 GInt / s (c) MI 250 3236439631283 50 100 150 200 𝑁 GInt / s (d) GPU Max 1100 Figure 2: Interactions per second vs number of particles Non the investigated architectures by the different implementations. Note the different scale of the y-axis for GH200 (b). We do not include the results for OpenACC on the MI 250 since they are about a factor of 100 lower than those from the other programming models and nearly independent of system size. might show better performance, when using algorithms more dependent on complex memory access patterns, e.g., short-range interactions, as in the case for Ewald Sum and the Particle-Particle Particle-Mesh algorithms. MultiXscale Deliverable 2.5 Page 12 numbers of nodes used, while for the larger number of nodes used the r-space runtime becomes dominant. Overall a parallel efficiency of ≈35% was reached on 128 GPUs (32 nodes). As a consequence of the finding a larger allocation of resources for the k-space part might be obvious to improve the total run-time. Therefore, the benchmark presented in Figure 5used three of four GPUs per node for the computation of the k-space contribution and one GPU of each node for the r-space contribution. However, it can be seen that the overall parallel efficiency decays faster than in the first benchmark. This is due to the fact that the r-space computation lost performance and dominates the overall parallel efficiency. The improved scaling behavior of the k-space computation can not be utilized, as the r-space computation always takes longer to execute. In a third benchmark a larger system was used to investigate the impact of system size on strong scalability. The results are shown in Figure 6. In this benchmark the resources were evenly split between both parts, rand k-space contributions. It can be seen that an efficiency of ≈60% was achieved. For all benchmarks the different computation parts are executed on sub-communicators which do not share processes. This means that they run independently in parallel but with different individual parallel efficiencies. As a result the sum of both parts is larger than the total runtime. First analysis indicates that this approach of splitting compute resources can be beneficial for node counts where the overall parallel efficiency starts to degrade. In order to be able to better guide the community in that respect further in-depth analysis will be performed. Figure 4: Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 16 million randomly placed particles in a periodic system of size 1003on Juwels Booster (Forschungszentrum Jülich). The resources are evenly split between the r-space and k-space computation. MultiXscale Deliverable 2.5 Page 13 Figure 5: Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 16 million randomly placed particles in a periodic system of size 1003on Juwels Booster (Forschungszentrum Jülich). Three of four GPUs are employed for the computation of the k-space contribution, the other GPUs are used for the r-space contribution. Figure 6: Strong scaling benchmark of the Ewald summation solver implemented in the library using a system of 32 million randomly placed particles in a periodic system of size 1263on Juwels Booster (Forschungszentrum Jülich). The resources are evenly split between the r-space and k-space computation. MultiXscale Deliverable 2.5 Page 14 5 Integration of the ALL load balancing library into LAMMPS The A Load balancing Library (ALL) is a modular, open-source software developed by the Jülich Supercomputing Centre to address the critical challenge of load imbalance in domain-decomposed simulations. In high-performance computing (HPC), particularly in particle-based methods like molecular dynamics or lattice-Boltzmann simulations, spatial domain decomposition is a common strategy to distribute computational workload. However, as simulations evolve, non-uniform distributions of particles or varying local computational costs can lead to significant load imbalances, reducing overall efficiency. ALL provides a suite of dynamic load-balancing strategies specifically designed to redistribute work across processors in such simulations, ensuring better utilization of HPC resources. By supporting multiple balancing schemes such as histogram and iterative grid partitioning ALL offers flexible and scalable solutions tailored to different application needs. Its integration into widely-used simulation codes like HemeLB and DL_MESO_DPD highlights its effectiveness in improving performance on large-scale systems, including GPU accelerated platforms. Within the scope of the current Deliverable, we highlight the most efficient and used methods, which are present in the library or being already implemented in LAMMPS: Tensor Decomposition (Local) This method splits the simulation domain along Cartesian axes using a fixed number of partitions in each direction (e.g., a grid layout). It is an iterative method that, for every dimension, reduces the work over the other dimensions and shifts the borders between two neighbours based on their cumulative work, making it simple but less adaptive to highly uneven load distributions. •Classic For the classic approach the work per dimension is collected using the average over the domains in that dimension. •Max The max variant considers the maximum work in each of the slices in cartesian direction and considers as objective function an even distribution of the maxima in the slices of each dimension. Staggered Grid The staggered grid is a special domain decomposition scheme that recursively cuts in each dimension in a defined order. It is a special case of the tiled decomposition as used for recursive bisection and is, in principle, able to produce a perfect decomposition of even load in the system. •Global (Global + Histogram) This method divides the system across domains into grid regions and adjusts the boundaries using global workload histograms. It optimizes the load balance by analyzing the distribution of computational work across the entire domain, then shifts domain boundaries accordingly. It is a global method and uses histogram-based refinement to redistribute work more evenly. •Local (Local + Iterative) Unlike the global variant, this method balances workloads based on local neighbor-toneighbor interactions. It iteratively adjusts boundaries by exchanging workload data with neighboring domains, gradually improving balance through a series of local updates. This local method uses iteration-based refinement, making it scalable and responsive to localized imbalances. Recursive Bisection (Global + Iterative) This is a global and hierarchical method that recursively splits the domain into two subregions at each step, aiming to evenly distribute the total workload. The splitting process can be guided by either workload histograms or iterative evaluation of load imbalance. It is effective for complex or non-uniform load distributions but can introduce communication overhead if applied too frequently. 5.1 ALL interface with LAMMPS ALL balancing methods have been integrated into LAMMPS using the standard definitions of so-called fixes. Once LAMMPS is build with -D PKG_ALLPKG=yes, the following options become available: fix ID group-ID balance/all keyword args ... required keyword/arg pairs every arg = nevery nevery = perform dynamic load balancing every this many steps grid style args = define grid style = staggered or tensor staggered args = none tensor args = max or classic max = use TENSOR_MAX method of ALL classic = use TENSOR method of ALL MultiXscale Deliverable 2.5 Page 15 weight style args = use weighted particle counts for the balancing style = group or neigh or time or var or store group args = Ngroup group1 weight1 group2 weight2 ... Ngroup = number of groups with assigned weights group1, group2, ... = group IDs weight1, weight2, ... = corresponding weight factors neigh factor = compute weight based on number of neighbors factor = scaling factor (> 0) time factor = compute weight based on time spend computing factor = scaling factor (> 0) var name = take weight from atom-style variable name = name of the atom-style variable store name = store weight in custom atom property defined by fix property/atom command name = atom property name (without d_ prefix) At least one of the following keyword/arg pairs is required. Only one (or none) load balancer is used in a load balancing step. global args = trigger_threshold bin_width trigger_threshold = float float = apply load balancer if the imbalance of the load is above this threshold bin_width = approximate width of bin (arbitrary positive value) local args = trigger_threshold (use local balancing) trigger_threshold = float float = apply load balancer if the imbalance of the load is above this threshold optional keyword arg pairs verbose args = none (verbose output to log/screen) file arg = filename filename = write stats per rank to filename In addition to the fix balance/all method, a specialized communication style was added in order to correctly and efficiently manage the MPI communicators set by the domain decomposition. Some examples of the available balancing methods are: # histogram method (staggered global), weighted by number of atoms, # with bin width 0.003, when imbalance >= 1.05 comm_style staggered zyx fix 5 all balance/all every 1000 grid staggered global 1.05 0.003 # histogram method (staggered global), weighted by time, with bin width 0.003, when imbalance >= 1.05 comm_style staggered zyx fix 5 all balance/all every 1000 grid staggered weight time 1.0 global 1.05 0.003 # staggered local method, weighted by number of atoms, when imbalance >= 1.05 comm_style staggered zyx fix 5 all balance/all every 100 grid staggered local 1.05 # tensor method, minimizing average work, weighted by number of atoms, when imbalance >= 1.05 comm_style brick fix 5 all balance/all every 50 grid tensor classic local 1.05 # tensor-max method, minimizing maximal work, weighted by number of atoms, when imbalance >= 1.05 comm_style brick fix 5 all balance/all every 50 grid tensor max local 1.05 MultiXscale Deliverable 2.5 Page 16 5.2 General LAMMPS Benchmarks With the intention to test LAMMPS performance on the most recent architecture present at the EuroHPC latest supercomputer JUPITER, we have run several benchmarks of LAMMPS with standard systems. The test consisted in weak and strong scaling on JEDI and JUWELS Booster systems both available at Forschungszentrum Jülich2. For the validation benchmarks, LAMMPS was compiled using GCC 12.3.0, CUDA 12.2, and OpenMPI 4.1.6, and the chosen FFT backend was cuFFT. In Figures 7and 8we present the results from porting LAMMPS to the newer architectures. (a) LJ (b) Rhodopsin Figure 7: Single node weak scaling of LAMMPS using JEDI and JUWELS Booster (Forschungzentrum Jülich) for the LJ and Rhodopsin benchmarks in LAMMPS using Kokkos. Scaled LJ simulation containing 55M atoms, and Scaled Rhodopsin simulation containing 16M atoms both with 1 MPI process per GPU. No results are presented for Juwels Booster with 1 and 2 GPUs for the Rhodopsin benchmark due to memory limitation to accommodate the data. As it can be observed in Figure 7both systems and architecture combinations show close to ideal scaling within single node configurations. The expected speed up in the newer architecture (2x) is achieved and the larger memory allows for the simulation of bigger system sizes with less computing resources. In Figure 7b one can observe that due to memory and communication characteristics, it is possible to achieve similar results for the Rhodopsin benchmark using one H100 of the GH200 node instead of three A100 GPUs. (a) LJ (b) Rhodopsin Figure 8: Strong scaling of LAMMPS using JEDI for the LJ and Rhodopsin benchmarks in LAMMPS using Kokkos. Strong scaling results present in the grouped bars of Figure 8show that for the LJ benchmark the minimum system size required to profit from allocating multi-node jobs was 55.4 million particles. Once taking into consideration the more complex system, i.e., the Rhodopsin benchmark, one may expect 16.5 million particles to be sufficient at the tested range, which seems to result from the more complex interactions. In fact, the Rhodopsin benchmark includes all electrostatic forces, bonds, dihedrals and improper contributions which results in a much larger computational work compared to the simple LJ case. 2The JEDI system has now been decommissioned and has been integrated into the JUPITER system. MultiXscale Deliverable 2.5 Page 17 Figure 9: Combined strong and weak scaling of LAMMPS at the pre-exascale supercomputer JEDI (Forschungszentrum Jülich). Solid bars use the straggered global load balancing method from the ALL library integrated into LAMMPS and hatched bars use the recursive bisection load balancing method native to LAMMPS. 5.3 Benchmark with balancing methods For testing the performance of the load balancing library, a synthetic benchmark was generated removing atoms from a spherical region at the centre of the LJ standard benchmark. The system was scaled in the same manner as the ones presented in the previous section. In Figure 9results show comparisons between load balancing methods, i.e., RCB, which is one of LAMMPS’s native benchmark methods, and the staggered global methods (ALL) using the same interval for checking the load imbalance and activating the load balancing procedure with a common criterion. The balancing commands added to the input script were: comm_style staggered zyx fix 5 all balance/all every 50 grid staggered global 0.9 0.003 in the case of using the method available through the integration with ALL, and : comm_style tiled fix 5 all balance 50 rcb 0.9 for LAMMPS’s recursive bisection algorithm (RCB). It can be observed that for all numbers of nodes the staggered/global method (solid bars) resulted in a performance improvement of 30% to 50% compared to RCB (hatched bars), which shows a promising capability of enhancing the performance of load imbalanced simulations on exascale supercomputers. These simulations were performed using the GPU package in LAMMPS through the use of -sf gpu command line. MultiXscale Deliverable 2.5 Page 18 6 Conclusion and Perspectives This deliverable presents the developments achieved through the efforts of Task 2.3. Sections 2to 5cover performance portability, electrostatic methods and implementation considerations, library API and benchmarks, and the integration of the ALL load balancing library into LAMMPS. These individual topics are interconnected, with portability aspects influencing the design choices of electrostatic solvers and their impact on target systems for optimal load balancing methods that interface with particle simulations. Kokkos has proven to be a reliable choice for porting codes to newer architectures, offering significant advantages in complex data structures and widespread community adoption. Although it may not yield the best overall P ¯ P ¯values, it consistently performs well on systems with a robust backend software stack. As demonstrated in Section 2, Kokkos matches the performance of native/vendor-recommended platform-dependent solutions without requiring laborious compiler optimizations. However, it is essential to note that platform-related compiler optimizations may positively impact Kokkos’ performance, and system maintainers are encouraged to explore this topic. Section 3presents algorithmic considerations aimed at extracting the desired performance from modern exascale architectures. Each independent step can utilize specific sub-communicators and run in parallel on CPU+GPU combinations or across modular partitions. Both Ewald summation and P3M methods employ either particle or mesh-based parallelism through Cabana methods or Kokkos Multi-Dimensional Range loops. The tentative API of the library and preliminary benchmarks are presented in Section 4. These benchmarks validate and highlight the library’s modular capability, demonstrating the splitting of rand k-space methods and their respective scalabilities for the Ewald summation method (Figures 4to 6). These benchmarks were performed by using two independent allocations for the rand k-space parts on Juwels Booster (Forschungszentrum Jülich). The discussion includes a critical assessment of the requirements and capabilities of a balancing method to run rand k-space parts most efficiently. Strong scaling of 16 and 32 million particles scaled from 4 to 128 GPUs reveals that parallel efficiency is affected by the manner in which computational effort of rand k-space parts is distributed onto resources. Initial analysis suggests that this approach to resource allocation can be beneficial for node counts where overall parallel efficiency begins to degrade. Further in-depth analysis will be conducted to provide more comprehensive guidance. The P3M implementation is currently undergoing validation, and its extraction to the new electrostatic library has commenced. We plan to analyze the impact of modular allocation with mesh methods using the same methodology presented for the Ewald Summation. Both methods will be optimized for running on Europe’s largest supercomputer, JUPITER. Lastly, the LAMMPS porting to the newest architecture has exhibited expected weak and strong scaling behavior. Performance on GH200 nodes is expected to be twice as fast as on the previous Nvidia GPU architecture (A100). For very large simulations, energy efficiency improves more than twofold, as fewer nodes and GPUs are required to perform equally sized simulations. Following the analysis of LAMMPS scaling behavior on both JEDI and Juwels Booster, preliminary results using the integration of ALL into LAMMPS are presented in Section 5.3. As the performance of load balancing depends on initial configuration and computing costs (described in Section 5.2), a more comprehensive study will be conducted to evaluate the differences between LAMMPS’ implemented balance methods and those added by the integration with the ALL library. The standardization of the interface is near to completion and will be made available to the community by creating a pull request to LAMMPS’ main repository. MultiXscale Deliverable 2.5 Page 19 Acknowledgements We thank the Jülich Supercomputing Centre, Forschungszentrum Jülich for access to the JURECA-DC Evaluation Platform, JUWELS Booster and JEDI. This work has received funding from the European High Performance Computing Joint Undertaking under the grant agreement 101093169 and was supported by the German Federal Ministry of Education and Research (BMBF) and the Ministry of Culture and Science (MKW) of the state of North-Rhine-Westphalia through funding of the Gauss Centre for Supercomputing (GCS). References Acronyms used HPC High Performance Computing FFT fast Fourier transform EuroHPC European High Performance Computing Joint Undertaking GPU Graphical Processing Unit MD Molecular Dynamics GCS Gauss Supercomputing Centre API Application Program Interface Software mentioned CUDA Compute Unified Device Architecture heFFTe highly efficient Fast Fourier Transform for exascale ALL A Load balancing Library ScaFaCoS Scalable Fast Coulomb Solvers URLs referenced Page ii https://www.multixscale.eu ... https://www.multixscale.eu https://www.multixscale.eu/deliverables ... https://www.multixscale.eu/deliverables Internal Project Management Link ... https://github.com/multixscale/planning/ [email protected] ... mailto:[email protected] http://creativecommons.org/licenses/by/4.0 ... http://creativecommons.org/licenses/by/4.0 Page 1 TOP500 ... https://www.top500.org Page 7 ScaFaCoS ... http://www.scafacos.de/ Page 9 Cabana ... https://github.com/ECP-copa/Cabana #799 ... https://github.com/ECP-copa/Cabana/pull/799 #800 ... https://github.com/ECP-copa/Cabana/pull/800 Page 18 JUPITER ... https://www.top500.org/lists/top500/2025/06/ Citations [1] T. Deakin, S. McIntosh-Smith, J. Price, A. Poenaru, P. Atkinson, C. Popa, and J. Salmon, “Performance portability across diverse computer architectures,” in 2019 IEEE/ACM International Workshop on Performance, Portability and Productivity in HPC (P3HPC). IEEE, 2019, pp. 1–13. [2] J. H. Davis, P. Sivaraman, I. Minn, K. Parasyris, H. Menon, G. Georgakoudis, and A. Bhatele, “Taking gpu programming models to task for performance portability,” arXiv preprint arXiv:2402.08950, 2024. [3] C. Phuong, N. Saied, and C. Tanis, “Assessing Kokkos performance on selected architectures,” in Latin American High Performance Computing Conference. Springer, 2019, pp. 170–184. [4] A. S. Dufek, R. Gayatri, N. Mehta, D. Doerfler, B. Cook, Y. Ghadar, and C. DeTar, “Case study of using Kokkos and SYCL as performance-portable frameworks for Milc-Dslash benchmark on NVIDIA, AMD and Intel GPUs,” MultiXscale Deliverable 2.5 Page 20 in 2021 International Workshop on Performance, Portability and Productivity in HPC (P3HPC), Nov. 2021, pp. 57–67. [5] M. Breyer, A. Van Craen, and D. Pflüger, “A comparison of SYCL, OpenCL, CUDA, and OpenMP for massively parallel support vector machine classification on multi-vendor hardware,” in Proceedings of the 10th International Workshop on OpenCL, ser. IWOCL ’22. New York, NY, USA: Association for Computing Machinery, May 2022, pp. 1–12. [6] Y. Ding, C. Xu, H. Qiu, Q. Wang, W. Dai, Y. Lin, and Y. Che, “Evaluating performance portability of SYCL and Kokkos: A case study on LBM simulations,” in 2023 IEEE Intl Conf on Parallel & Distributed Processing with Applications, Big Data & Cloud Computing, Sustainable Computing & Communications, Social Computing & Networking (ISPA/BDCloud/SocialCom/SustainCom), Dec. 2023, pp. 328–335. [7] W.-C. Lin, T. Deakin, and S. McIntosh-Smith, “Evaluating ISO C++ parallel algorithms on heterogeneous HPC systems,” in 2022 IEEE/ACM International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS). IEEE, 2022, pp. 36–47. [8] R. A. Bartolomeu, R. Halver, J. H. Meinke, and G. Sutmann, “Effect of implementations of the n-body problem on the performance and portability across gpu vendors,” Future Generation Computer Systems, vol. 169, p. 107802, 2025. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167739X25000974 [9] K. Asanovic, R. Bodik, B. C. Catanzaro, J. J. Gebis, P. Husbands, K. Keutzer, D. A. Patterson, W. L. Plishker, J. Shalf, S. W. Williams, and K. A. Yelick, “The landscape of parallel computing research: A view from Berkeley,” Electrical Engineering and Computer Sciences, University of California at Berkeley, Technical Report No. UCB/EECS-2006183, December, vol. 18, no. 2006-183, p. 19, 2006. [10] C. R. Trott, D. Lebrun-Grandié, D. Arndt, J. Ciesko, V. Dang, N. Ellingwood, R. Gayatri, E. Harvey, D. S. Hollman, D. Ibanez, N. Liber, J. Madsen, J. Miles, D. Poliakoff, A. Powell, S. Rajamanickam, M. Simberg, D. Sunderland, B. Turcksin, and J. Wilke, “Kokkos 3: Programming model extensions for the exascale era,” IEEE Transactions on Parallel and Distributed Systems, vol. 33, no. 4, pp. 805–817, 2022. [11] J. H. Davis, P. Sivaraman, I. Minn, K. Parasyris, H. Menon, G. Georgakoudis, and A. Bhatele, “An Evaluative Comparison of Performance Portability across GPU Programming Models,” arXiv preprint arXiv:2402.08950, 2024. [12] W.-C. Lin, S. McIntosh-Smith, and T. Deakin, “Preliminary report: Initial evaluation of StdPar implementations on AMD GPUs for HPC,” Jan. 2024. [13] A. Arnold, F. Fahrenberger, C. Holm, O. Lenz, M. Bolten, H. Dachsel, R. Halver, I. Kabadshow, F. Gähler, F. Heber, J. Iseringhausen, M. Hofmann, M. Pippig, D. Potts, and G. Sutmann, “Comparison of scalable fast methods for long-range interactions,” Phys. Rev. E, vol. 88, p. 063308, Dec 2013. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.88.063308 [14] M. Deserno and C. Holm, “How to mesh up ewald sums. i. a theoretical and numerical comparison of various particle mesh routines,” The Journal of Chemical Physics, vol. 109, no. 18, pp. 7678–7693, 11 1998. [Online]. Available: https://doi.org/10.1063/1.477414 [15] R. W. Hockney and J. W. Eastwood, Computer simulation using particles. crc Press, 2021. [16] J. J. Cerdà, V. Ballenegger, O. Lenz, and C. Holm, “P3m algorithm for dipolar interactions,” The Journal of Chemical Physics, vol. 129, no. 23, p. 234104, 12 2008. [Online]. Available: https://doi.org/10.1063/1.3000389 [17] S. Slattery, S. T. Reeve, C. Junghans, D. Lebrun-Grandié, R. Bird, G. Chen, S. Fogerty, Y. Qiu, S. Schulz, A. Scheinberg, A. Isner, K. Chong, S. Moore, T. Germann, J. Belak, and S. Mniszewski, “Cabana: A performance portable library for particle-based simulations,” Journal of Open Source Software, vol. 7, no. 72, p. 4115, 2022. [Online]. Available: https://doi.org/10.21105/joss.04115