# SOMA-BD: Brownian dynamics simulation for soft matter on GPU

Benoît Crespin, Heinich Porro, Manuella Cerbelaud, Arnaud Videcoq, Julien Gerhards

► To cite this version:

Benoît Crespin, Heinich Porro, Manuella Cerbelaud, Arnaud Videcoq, Julien Gerhards. SOMA-BD: Brownian dynamics simulation for soft matter on GPU. Engineering with Computers, 2024, 10.1007/s00366-024-02072-1. [hal-04740841](https://hal.science/hal-04740841v1)

HAL Id: hal-04740841

https://hal.science/hal-04740841v1

Submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SOMA-BD: Brownian Dynamics Simulation for
Soft Matter on GPU

Benoit Crespin¹*, Heinich Porro¹,², Manuella Cerbelaud²,
Arnaud Videcoq², Julien Gerhards²

¹Univ. Limoges, CNRS, XLIM, UMR 7252, F-87000, Limoges, France.

²Univ. Limoges, CNRS, IRCER, UMR 7315, F-87000, Limoges, France.

*Corresponding author(s). E-mail(s): benoit.crespin@unilim.fr;

**Abstract**

This paper introduces a high-performance Brownian Dynamics simulation code designed for Nvidia GPUs. The code prioritizes accuracy with double-precision calculations and utilizes a bucket sort neighbor list for efficient particle inter-action management. Additionally, a CUDA occupancy-aware autotuning system optimizes thread block size for diverse hardware configurations. Benchmarks demonstrate significant performance gains compared to existing tools, partic-ularly for large-scale simulations with various interaction potentials (Lennard-Jones, Yukawa, DLVO). The software offers user-friendly configuration via JSON files and supports restarting simulations from both binary and XYZ data for-mats. Future work includes expanding supported interaction potentials, refining autotuning and investigating scalability. This open-source software empowers researchers in material science and related fields to conduct efficient and accurate BD simulations.

**Keywords:** GPU Computing, Brownian Dynamics, Colloidal Suspensions

001

002

003

004

005

006

007

008

009

010

011

012

013

014

015

016

017

018

019

020

021

022

023

024

025

026

027

028

029

030

031

032

033

034

035

036

037

038

039

040

041

042

043

044

045

046
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 1 Simulation of a colloidal system with two types of particles of varying radii with a Yukawa potential (image rendered with Ovito Stukowski (2010)).</figcaption></figure>

# 1 Introduction

The combination between particle interactions and thermal fluctuations gives rise to distinctive features of soft matter, which includes a wide range of materials such as liquids, colloids, and polymers. In materials science and related subjects, it is essential to comprehend their behavior and computer simulations are frequently required due to the complexity of real-world soft matter systems. Simulations of Brownian dynamics (BD) describe the motion of small particles subjected to Brownian motion, such as colloids immersed in a liquid. Collisions with the surrounding fluid molecules, which are too many to precisely simulate, cause this erratic, jittery movement. For soft matter research, BD offers a considerable efficiency advantage over standard molecular
dynamics (MD) simulations. This is because it treats the solvent implicitly, concentrating only on the interactions and dynamics of the colloids, which makes it especially useful for modeling big and intricate soft matter systems. Figure 1 illustrates the aggregation process most often observed, showing a typical image generated after simulating a colloidal system with two types of particles of varying radii.

Specialized processors called graphics processing units, or GPUs, were first created to speed up the rendering of visuals. Their parallel processing design, which enables them to handle several tasks at once, is actually their greatest strength. For computationally demanding scientific simulations like BD, this makes them perfect. Although GPUs have enormous performance potential, there are obstacles to be addressed. To guarantee precise and effective BD simulations on the GPU, certain techniques are needed for example to manage neighbor lists — which monitor interacting particles — efficiently.

This paper presents a novel GPU-based software specifically designed for performing BD simulations of colloidal suspensions. Our program tackles a number of important issues when employing BD on GPUs to simulate colloidal suspensions. First, it addresses how to use GPU resources effectively for simulations with a modest number of colloids. This is made possible by certain optimizations that guarantee the GPU is always in use. Second, the software offers flexibility in particle size representation. In order to accommodate situations with heterogeneous colloidal systems, it enables researchers to establish a size distribution or specify individual particle radii. Moreover, simplicity and maintainability are prioritized in the C++ writing of the code. This design decision enables material scientists, who may not be proficient in GPU programming, to easily adapt and expand the code to suit their particular research requirements. Last but not least, the program maintains double precision for particle position storage, guaranteeing excellent simulation accuracy and precision — a need for a number of material science applications.
The paper is organized as follows. In section 2, we present an overview of the lit-
erature on BD simulations that has already been published, emphasizing the present
state of the art in terms of GPU-based software availability. The general organiza-
tion of our software is then covered in section 3, including how the code operates
and which GPU optimizations are applied to make sure that GPU resources are used
as efficiently as possible. Section 4 includes the results that were achieved using our
software, together with a detailed commentary that highlights the features and pos-
sible applications. In order to enable researchers to customize the software to their
own requirements, we describe in section 5 the implementation specifics, such as the
important design decisions that encourage code simplicity.

## 2 Existing works

Molecular Dynamics (MD) software packages may be used to conduct Brownian
Dynamics simulations. Numerous systems are available in computational physics
and chemistry, such as HOOMD Anderson et al (2020); Howard et al (2016, 2019),
LAMMPS Thompson et al (2022), GROMACS Abraham et al (2015); Kim et al (2023),
TinkerHP Adjoua et al (2021a,b), HALMD Colberg and Höfling (2011) or UAMMD
Hashemi et al (2023). By utilizing these broad MD frameworks’ flexibility and effi-
ciency, scientists may easily simulate the stochastic motion of particles immersed in a
fluid as a result of thermal fluctuations, simulating the behavior of biological molecules
in solution, polymers, and colloidal particles. Scientists can now study complex sys-
tems at the molecular scale with unprecedented accuracy and detail thanks to the
integration of Brownian Dynamics capabilities into these well-established simulation
environments.

All of the aforementioned software programs can operate on GPUs, however con-
figuring them for basic BD simulations can occasionally be difficult. Compiling the
software with the necessary libraries and dependencies is usually the first step in
configuring GPU support, and it may call for some technical know-how. The setup pro-
cedure may also get more complicated if simulation settings are optimized for effective
GPU use and compatibility with certain hardware setups. Despite these difficulties,
researchers looking to improve the speed and effectiveness of their BD simulations
should still pursue GPU acceleration due to its performance advantages, which include
quicker simulation times and greater scalability for bigger systems. In addition, users
may find tutorials, documentation materials, and community support forums to help
them with the installation, compilation, and usage of these software packages on
GPU-accelerated platforms.

Neighbor lists, also called Verlet lists, are essential for effectively computing particle
interactions in different types of particle-based simulations, including BD Winkler et al
(2018); Lipscomb et al (2012); Tran et al (2018). These lists are crucial for determining
particle pairings that are inside a specific interaction range. This helps to avoid the
computationally costly effort of analyzing interactions between every pair of particles.
Neighbor lists are essentially designed to be updated on a regular basis depending on
particle motions during the simulation. At first, any particle that is within a given
cutoff distance is added to the list and regarded as a neighbor. Particles in the system
move as the simulation goes on, possibly changing how close they are to one another. If
these changes are too important, the neighbor lists are updated to reflect these changes
in particle locations and preserve efficiency and accuracy. This makes BD simulations
computationally practical for researching huge systems over long extended timescales.

Effectively updating Verlet lists is essential to accomplishing high-performance
simulations. One popular strategy is the grid-based technique, also referred to as
cell-lists Hoetzlein (2014); Gross et al (2019); Chisholm et al (2020); Bender, Jan
(2024). Using this approach, the simulation domain is partitioned into a regular grid
of cells, and the locations of the particles are used to assign them to the cells. Particles
can drastically reduce the processing cost of identifying possible interaction partners
by navigating adjacent cells within a predetermined radius. Hierarchical techniques
Fernández-Fernández et al (2022); Howard et al (2016, 2019); Ohno et al (2017) use
effective tree traversal algorithms to arrange particles into hierarchical data structures
that facilitate quick neighbor searches. Because hierarchical methods are flexible and
scalable, they may be used in simulations with different system densities and sizes.
Finally, recent developments take advantage of hardware-accelerated ray tracing to
effectively find particle neighbors by speeding up closest neighbor searches using real-
time ray tracing capabilities Zhu (2022). These techniques yield impressive advances
in neighbor list creation performance but require the newest characteristics of GPU
architectures such as NVIDIA's RT Cores. In GPU-accelerated simulations, the devel-
opment of effective neighbor list creation methods is still an important research topic,
however grid-based strategies are usually thought to be the most practical.

Managing the size of neighbor lists in memory is also essential for optimizing
simulation performance and resource utilization. If the allocated memory for neighbor
lists is too large, it can lead to unnecessary memory consumption and inefficient use of
resources. Conversely, if the memory allocation is too small, it may not accommodate
the actual number of effective neighbors for each particle, potentially resulting in
inaccurate simulations Kim et al (2023). Setting this option can be challenging, for
example the number of neighbors is restricted to 70 in Bender, Jan (2024) but can be
modified in the code.

Other possible optimizations include auto-tuning Hijma et al (2023), defined as
a way to optimize the GPU thread block sizes for each kernel according to the
GPU architecture that the code is running on. With this method, thread block sizes
are dynamically adjusted to enhance performance across various GPU architectures.
Additionally, other research focuses on memory cache optimizations, although they
typically apply to scenarios involving small numbers of particles Dice et al (2018).
# 3 Code operation and GPU optimizations

277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322

Figure 2 summarizes the workflow of our simulation method. Particle positions are either imported from an input file or defined randomly during the setup phase. Simulation parameters include the number of iterations, the number of particles and their density, their size distribution and the type of interacting forces. The list of forces available and a description of the restart procedure mentioned on Fig. 2 are provided in section 5.

Our solution uses neighbor lists in the main simulation loop, as other existing software presented previously. These lists are updated when particle displacements exceed a given threshold. The updating process is based on bucket sorting, a grid-based method similar to Bender, Jan (2024) and adapted to periodic conditions. At each iteration step, the neighbor lists are used to determine the interacting particles, compute their interaction forces and calculate their displacement at the next iteration. Every one of these tasks relates to a particular GPU kernel, which are optimized for maximum efficiency prior to the simulation. The subsequent sections provide an explanation of this optimization procedure as well as specifics on the neighbor list maintenance and bucket sorting techniques.

## 3.1 Bucket sort

Our code implements a GPU-based bucket sort algorithm, following the same grid-based approach described in Gross et al (2019); Hoetzlein (2014); Bender, Jan (2024) and specifically designed for fixed radius neighbor search with periodic conditions.

At the beginning of the algorithm, the grid structure is initialized on the GPU. The algorithm proceeds by inserting each particle into the appropriate cell in the grid. Once all particles have been assigned to one cell, a sorting operation is performed based on the cell index. The entire algorithm, from particle insertion to sorting and neighbor search, is designed to take advantage of the parallel processing capabilities of
<figure><img src="image_4.png" /><figcaption>Fig. 2 Overview of our simulation workflow</figcaption></figure>

the GPU, reduce memory latency and optimize cache utilization, resulting in improved overall performance.

The approach of storing the particle count and the index of the first particle in each cell is a commonly used and efficient method in most grid-based implementations on the GPU. Our implementation relies on a slightly different idea, using a single `int2` array representing a pair of values for each cell: the index of the first particle in the cell, and the index of the last one. This array can be filled with a single GPU kernel iterating over the sorted particle array in parallel, using atomic min and max operations to update the first and last indices for the corresponding cell in the array (see Figure 3). In addition to minimizing memory usage, the use of a single buffer
<figure><img src="image_5.png" /><figcaption>**Fig. 3** GPU arrays built from an initial configuration of particles and grid cells (left): arrays storing particle counts and indices of the first particle in each cell (middle) are replaced by a single array storing the indices of the first and last particles (right). For the sake of clarity, the Z-order curve is not used to assign hash values to cells in this figure.</figcaption></figure>

enables contiguous storage of data for a given cell, encompassing both the index of the initial particle and the index of the last particle. This contiguous storage layout improves cache utilization because when accessing neighboring particles, the data for a specific cell is loaded into the cache as a single block. Our experiments showed that this faster retrieval of neighboring particles and more efficient processing of neighbor search operations can improve performance by around 10%.

Bucket sort does have a few drawbacks such as imbalanced workloads across the grid cells, particularly when the particle distribution is not uniform, which poses a problem because the traditional bucket sort assumes a one-to-one mapping between particles and grid cells. Alternative approaches such as hierarchical structures Fernández-Fernández et al (2022) may offer more optimized solutions if the differences between the various interaction radii are very large. Our solution is simply to use the largest radius size as the grid cell size. In this approach, each grid cell can fully accommodate the largest particle and avoids incomplete particle coverage within the grid cells, as no particle spans several cells. This may result in a slight lack of efficiency and increased memory usage, but remains largely easier to implement than a hierarchical approach.

369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
## 3.2 Managing neighbour lists

The Verlet list of each particle is constructed by identifying their neighbors within a cutoff distance $R_l$, which is typically obtained by multiplying the interaction radius between particles $R_c$ by an increasing factor. This factor called "Rc_multiplicator" is a parameter of the simulation and must be carefully chosen: if too large, then the amount of potential neighbours may result in a larger number of unnecessary distance computations at each iteration. On the other hand, if $R_l$ is too close to $R_c$ then the neighbour lists must be reconstructed more often, which again results in degraded performances although the simulation will still produce correct results. Several experiments are necessary to determine an optimal value, but our experiments show that a value of 1.1 usually results in acceptable performances.

In order to appropriately reflect possible interactions between particles, the Verlet lists must be updated when particles move and their locations change during the simulation. The decision criterion for recalculating neighbor lists is, in fact, a crucial component of Verlet list-based techniques. We employ the test described in Tran et al (2018), in which we determine if any particle has migrated farther than the threshold distance $(R_l - R_c)/2$ for the largest particle since the latest update. If that is the case, we use the bucket sort technique from the previous section to recalculate the complete set of neighbor lists.

Some authors have questioned whether the use of neighbor lists is really worthwhile, since a certain amount of memory is required to store these lists Winkler et al (2018). Indeed, an alternative strategy might be to use the bucket sort method at each iteration, calculating interactions on the fly without the need to store neighbors in memory. This could apply to cases where each particle has a small number of neighbors on average, or where the high speed of the particles means that the test described in the previous paragraph leads to lists being reconstructed almost every iteration. However, in the case of molecular dynamics simulations, the general case is that of low
speeds and a neighborhood that does not change much between two computational iterations. In most of our simulations, this means that the lists are recalculated only every 50 to 200 iterations, making the principle of neighbor lists much more beneficial than a simple grid.

Another improvement could be to exploit the symmetry of particle interactions, to avoid, for example, if a particle A is interacting with a particle B, storing B in A's neighbor list and also A in B's neighbor list. However, we considered that the main consequence of this idea would make the code and memory pattern of neighbor lists more complex, which may in turn degrade performances Winkler et al (2018).

### 3.3 Automatic calculation of kernel sizes

Autotuning kernels on GPUs involves dynamically adjusting the parameters and configurations of GPU kernels to achieve optimal performance for a given hardware architecture and workload. This process typically includes tuning thread block sizes, memory usage patterns, and other kernel launch parameters, considering factors such as the GPU's architecture, the nature of the computation, and the size of the data being processed. Our study focused on thread block sizes, by systematically exploring different sizes from 32 to 2048 during the first 100 iterations of the simulation and keeping the one with the best results for each kernel. The choice of the first 100 iterations is motivated by the fact that particle behavior is relatively stable over the entire duration of our simulations; for other cases, it would probably be necessary to test the best configurations at regular intervals during the simulation to adjust the parameters if necessary.

For purely coalescent kernels, autotuning should give results comparable to those obtained using the `cudaOccupancyMaxPotentialBufferSize` routine in CUDA which calculates the optimal block size that maximizes the overall occupancy and efficiency of the kernel execution. This is the case for example for the kernels computing
the cell index for each particle or calculating the interaction force between neighbouring particles. For more complex kernels, such as those used to construct GPU arrays for each cell as in Figure 3, autotuning provides much better results. While `cudaOccupancyMaxPotentialBufferSize` focuses primarily on maximizing the occupancy by determining the optimal block size, it does not account for the specific computational patterns and memory access behaviors. This is illustrated in the next section, showing the difference in performances between optimized and non-optimized block sizes on different GPU architectures.

# 4 Performance measurements and discussion

We evaluated the performance of our software across a range of configurations, considering several critical parameters. The benchmarks were conducted on different types of Nvidia GPUs, from the oldest P100 (Pascal architecture) to the most recent H100 (Hopper), showcasing the software's adaptability to varying hardware capabilities. Table 1 presents various characteristics of GPUs. The table is sorted by the last column, which shows the theoretical performance for 64-bit floating-point operations, measured in GFlops.

<table><caption>**Table 1** Different types of Nvidia GPUs used for our benchmark measurements</caption><thead><tr><th>GPU</th><th>Architecture</th><th>Release date</th><th>CUDA Cores</th><th>Memory (GB)</th><th>FP64 (GFLOPS)</th></tr></thead><tbody><tr><td>1080Ti</td><td>GeForce 10</td><td>2017</td><td>3584</td><td>11</td><td>354.4</td></tr><tr><td>3090</td><td>GeForce 30</td><td>2020</td><td>10496</td><td>24</td><td>556</td></tr><tr><td>A40</td><td>Tesla Ampere</td><td>2020</td><td>10752</td><td>48</td><td>584.6</td></tr><tr><td>P100</td><td>Tesla Pascal</td><td>2016</td><td>3584</td><td>16</td><td>4763</td></tr><tr><td>L40</td><td>Tesla Ada</td><td>2022</td><td>18176</td><td>48</td><td>1414</td></tr><tr><td>H100</td><td>Tesla Hopper</td><td>2023</td><td>14592</td><td>80</td><td>25610</td></tr></tbody></table>

All the simulations presented below were carried out over 1 million iterations. The number of particles tested ranged from 10 thousands to 1 million, with different types of interaction forces. One widely used model is the Lennard-Jones (LJ) potential, which
captures the competition between attractive and repulsive forces. The LJ potential 553 between two particles with identical radii $a_c$ can be expressed mathematically as: 554

$$U_{LJ}(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \qquad (1)$$

where $\sigma = 2a_c$, $\epsilon$ is the depth of the potential well, and $r$ is the interparticle 560 distance Israelachvili (2011). This model can be modified to account for particles with 561 different radii using the Lorentz-Berthelot mixing rules. The corresponding benchmark 562 is denoted below as LJM, with a normal distribution of particle radii. 563

The Yukawa potential offers a versatile alternative that can capture short-range 567 interactions with a decaying tail. Here, the Yukawa potential is defined as: 568

$$U_Y(r) = \epsilon q_i q_j \frac{2a_c}{r} \exp \left( -\frac{r-2a_c}{\lambda} \right) \qquad (2)$$

where $\epsilon$ defines the well depth, $q_i$ the sign of the surface charge of particle $i$, $\lambda$ is the screening length controlling the decay of the interaction, and $r$ is the interparticle 574 distance Cerbelaud et al (2017). The Yukawa potential is not valid for a distance 575 shorter than the contact distance, for which a repulsive linear potential is used to 576 avoid interpenetration. 577

Finally, another common interaction potential is the Derjaguin-Landau-Verwey-582 Overbeek (DLVO) theory, which describes the combined effect of van der Waals 583 attraction and electrostatic repulsion between charged particles. The DLVO potential 584 depends on various factors like surface potential, ionic strength, and particle sizes, and 585 can be expressed as: 586

$$U_{DLVO}(r) = U_{vdW} + U_{electrostatic} \qquad (3)$$

with

598
599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

$$U_{vdW} = -\frac{A}{6} \left[ \frac{2a_i a_j}{r^2 - (a_i + a_j)^2} + \frac{2a_i a_j}{r^2 - (a_i - a_j)^2} + \ln \left( \frac{r^2 - (a_i + a_j)^2}{r^2 - (a_i - a_j)^2} \right) \right] \quad (4)$$

and

$$U_{electrostatic} = \pi \epsilon a_{eff} (\psi_i^2 + \psi_j^2) \left[ \frac{2\psi_i \psi_j}{\psi_i^2 + \psi_j^2} \ln \left( \frac{1 + \exp(-\kappa h)}{1 - \exp(-\kappa h)} \right) + \ln(1 - \exp(-2\kappa h)) \right] \quad (5)$$

where A is the Hamaker constant related to van der Waals attraction, a_i the radius of particle i, r the distance between particles, ϵ the dielectric constant of solvant, κ the inverse Debye length (κ = 1/λ), ψ_i the surface potential of the particle i, h the distance between the surface of particles i and j (h = r - a_i - a_j) and a_{eff} = a_i a_j / (a_i + a_j).

The DLVO potential becomes infinite at contact. To avoid divergence, this potential is cut at a potential well depth defined by the user and a repulsive linear potential is used for distances lower than the contact distance to avoid interpenetration.

Our measurements, summarized in Tables 2 and 3 (also shown on Figures 4 to 6 with a logarithmic scale), indicate that our software scales effectively with increasing particle numbers. It also scales well with an increase in density, representing the fraction of the simulation volume occupied by the volume of colloidal particles and set respectively at 5%, 10% and 20%. This parameter has a direct influence on the average number of neighbors that must be stored and processed for each particle during the simulation. Nevertheless, performance seems to be less density-dependent on the most recent architectures, as can be seen with the LJ simulation on the H100 card for 100K particles (see highlighted cells in Table 2). In this case, the density of 10% or 20% does not seem to have a particular influence, even if this behaviour is not found with a higher number of particles. This may be explained by the large amount of accessible memory and more optimized cache transfers.
<table><caption>**Table 2** Computation times in seconds for 1 million iterations with Lennard-Jones interactions (fixed size or normal distribution)</caption><thead><tr><th rowspan="2">GPU</th><th rowspan="2">Particles</th><th colspan="3">LJ</th><th colspan="3">LJM</th></tr><tr><th>5%</th><th>10%</th><th>20%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td rowspan="3">1080Ti</td><td>10K</td><td>74.09</td><td>96.36</td><td>103.48</td><td>69.18</td><td>84.55</td><td>95.66</td></tr><tr><td>100K</td><td>423.1</td><td>606</td><td>710.22</td><td>375.03</td><td>509.02</td><td>663.22</td></tr><tr><td>1M</td><td>3182.52</td><td>4283.91</td><td>5379.67</td><td>2626.74</td><td>3546.75</td><td>5058.53</td></tr><tr><td rowspan="3">A40</td><td>10K</td><td>59.98</td><td>69.02</td><td>74.13</td><td>51.4</td><td>61.34</td><td>72.04</td></tr><tr><td>100K</td><td>289.16</td><td>374.29</td><td>460.47</td><td>237.12</td><td>316.35</td><td>432.35</td></tr><tr><td>1M</td><td>2411.85</td><td>3207.02</td><td>4001.171</td><td>1978.19</td><td>2641.36</td><td>3744.98</td></tr><tr><td rowspan="3">3090</td><td>10K</td><td>51.82</td><td>62.97</td><td>70.88</td><td>49.6</td><td>55.38</td><td>67.49</td></tr><tr><td>100K</td><td>257.19</td><td>358.21</td><td>426.42</td><td>210.18</td><td>302.39</td><td>400.46</td></tr><tr><td>1M</td><td>2167.72</td><td>2943.97</td><td>3681.21</td><td>2230.05</td><td>2237.21</td><td>3444.42</td></tr><tr><td rowspan="3">P100</td><td>10K</td><td>70.88</td><td>78.25</td><td>77.85</td><td>66.63</td><td>70.98</td><td>71.08</td></tr><tr><td>100K</td><td>231.66</td><td>276.27</td><td>338.6</td><td>209.84</td><td>236.96</td><td>321.12</td></tr><tr><td>1M</td><td>1597.48</td><td>2058.87</td><td>2666.61</td><td>1411.97</td><td>1785.45</td><td>2586.11</td></tr><tr><td rowspan="3">L40</td><td>10K</td><td>39.4</td><td>52.76</td><td>52.22</td><td>41.74</td><td>47.16</td><td>55.36</td></tr><tr><td>100K</td><td>124.36</td><td>156.02</td><td>188.76</td><td>103.24</td><td>132.57</td><td>177.1</td></tr><tr><td>1M</td><td>1090.45</td><td>1432.69</td><td>1764.74</td><td>908.37</td><td>1183.67</td><td>1648.3</td></tr><tr><td rowspan="3">H100</td><td>10K</td><td>36.57</td><td>44.79</td><td>47.74</td><td>34.5</td><td>37.7</td><td>42.51</td></tr><tr><td>100K</td><td>50.93</td><td>59.78</td><td>60.6</td><td>45.89</td><td>51.38</td><td>58.91</td></tr><tr><td>1M</td><td>337.68</td><td>392.68</td><td>434.84</td><td>326.08</td><td>358.20</td><td>417.22</td></tr></tbody></table>

<table><caption>**Table 3** Computation times in seconds for 1 million iterations with Yukawa and DLVO interactions</caption><thead><tr><th rowspan="2">GPU</th><th rowspan="2">Particles</th><th colspan="3">Yukawa</th><th colspan="3">DLVO</th></tr><tr><th>5%</th><th>10%</th><th>20%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td rowspan="3">1080Ti</td><td>10K</td><td>92.87</td><td>109.32</td><td>120.77</td><td>141.38</td><td>183.04</td><td>202.84</td></tr><tr><td>100K</td><td>540.09</td><td>647.04</td><td>1010.49</td><td>811.4</td><td>1113.7</td><td>1452.67</td></tr><tr><td>1M</td><td>4233.68</td><td>5618.16</td><td>7521.85</td><td>6873.57</td><td>9630.1</td><td>12980.65</td></tr><tr><td rowspan="3">A40</td><td>10K</td><td>73.72</td><td>86.88</td><td>95.11</td><td>109.46</td><td>138.24</td><td>167.59</td></tr><tr><td>100K</td><td>369.36</td><td>478.29</td><td>604.91</td><td>593.3</td><td>817.72</td><td>1134.26</td></tr><tr><td>1M</td><td>3218.73</td><td>4234.54</td><td>5337.35</td><td>5178.78</td><td>7325.18</td><td>10290.75</td></tr><tr><td rowspan="3">3090</td><td>10K</td><td>60.91</td><td>73.35</td><td>83.66</td><td>89.65</td><td>113.52</td><td>144.32</td></tr><tr><td>100K</td><td>321.25</td><td>445.42</td><td>549.02</td><td>509.65</td><td>718.26</td><td>966.4</td></tr><tr><td>1M</td><td>2838.5</td><td>3730.59</td><td>4746.58</td><td>6309.13</td><td>6391.64</td><td>8826.74</td></tr><tr><td rowspan="3">P100</td><td>10K</td><td>77.49</td><td>71.04</td><td>80.08</td><td>85.48</td><td>93.65</td><td>112.46</td></tr><tr><td>100K</td><td>238.24</td><td>279.62</td><td>333.26</td><td>272.17</td><td>310.17</td><td>378.05</td></tr><tr><td>1M</td><td>1648.23</td><td>2009.08</td><td>2566.63</td><td>1840.07</td><td>2198.52</td><td>2688.54</td></tr><tr><td rowspan="3">L40</td><td>10K</td><td>59.48</td><td>65.35</td><td>69.76</td><td>77.88</td><td>89.08</td><td>115.19</td></tr><tr><td>100K</td><td>164.48</td><td>206.45</td><td>254.96</td><td>280.74</td><td>383.12</td><td>513.51</td></tr><tr><td>1M</td><td>1430.74</td><td>1870.87</td><td>2325.91</td><td>2293.15</td><td>3189.45</td><td>4386.64</td></tr><tr><td rowspan="3">H100</td><td>10K</td><td>38.24</td><td>41.61</td><td>45.42</td><td>46.44</td><td>52.28</td><td>67.33</td></tr><tr><td>100K</td><td>58.52</td><td>62.31</td><td>67.7</td><td>79.54</td><td>85.56</td><td>105.09</td></tr><tr><td>1M</td><td>373.47</td><td>417.71</td><td>464.79</td><td>474.75</td><td>570.06</td><td>726.41</td></tr></tbody></table>

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690
<figure><img src="image_6.png" /><figcaption>Fig. 4 Graphical representation of our benchmark measurements from Table 2 (LJ)</figcaption></figure>

Our software handles various interaction forces with high performance, and as expected our measurements show that calculation speed also depends on the complexity of the potential used. This is particularly true when the mathematical formulation involves the use of functions such as exp, which generally perform poorly on GPUs. One point to note in our case is the performance degradation on the L40 card, highlighted in table 3 for the DLVO potential. This behaviour could be due to the card's low capacity for double-precision calculation, which in this case would lead us to choose the P100 GPU even if it relies on an older architecture.

We conducted more tests to look at the effects of CUDA kernel size on computational performance for L40 and H100 GPUs. While the `cudaOccupancyMaxPotentialBlockSize` function consistently returned a value of 640 for both cards, our autotuning system described in section 3.3 identified significantly different optimal kernel sizes (shown on the first line of each experiment in Table 4). The second line of each entry in the table presents the acceleration factor achieved by utilizing the optimal kernel size compared to the non-optimal size.
<figure><img src="image_7.png" /><figcaption>Fig. 5 Graphical representation of our benchmark measurements from Table 2 (LJM)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6 Graphical representation of our benchmark measurements from Table 3 (Yukawa)</figcaption></figure>

As evident from Table 4, the acceleration factor is generally greater than 1, indicating faster computations with the optimal kernel sizes. However, there are some exceptions highlighted in the table that warrant further investigation. Understanding

737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
<figure><img src="image_9.png" /><figcaption>Fig. 7 Graphical representation of our benchmark measurements from Table 3 (DLVO)</figcaption></figure>

the reasons behind these specific cases where non-optimal sizes perform better will be crucial for further optimization.

<table><caption>Table 4 Optimal CUDA kernel sizes and acceleration obtained compared to non-optimal sizes on L40 and H100 cards</caption><thead><tr><th rowspan="2">GPU</th><th rowspan="2">Particles</th><th colspan="3">LJ</th><th colspan="3">LJM</th><th colspan="3">Yukawa</th><th colspan="3">DLVO</th></tr><tr><th>5%</th><th>10%</th><th>20%</th><th>5%</th><th>10%</th><th>20%</th><th>5%</th><th>10%</th><th>20%</th><th>5%</th><th>10%</th><th>20%</th></tr></thead><tbody><tr><td rowspan="4">L40</td><td>10K</td><td>96</td><td>32</td><td>32</td><td>32</td><td>32</td><td>32</td><td>32</td><td>96</td><td>32</td><td>96</td><td>96</td><td>32</td></tr><tr><td>100K</td><td>2.34</td><td>2.20</td><td>2.94</td><td>2.60</td><td>2.96</td><td>3.05</td><td>2.25</td><td>2.65</td><td>3.10</td><td>2.73</td><td>3.28</td><td>3.43</td></tr><tr><td>1M</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td></td><td></td><td>1.21</td><td>1.27</td><td>1.46</td><td>1.78</td><td>1.80</td><td>1.66</td><td>1.42</td><td>1.48</td><td>1.50</td><td>1.35</td><td>1.38</td><td>1.42</td></tr><tr><td rowspan="4">H100</td><td>10K</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>100K</td><td>1.21</td><td>1.27</td><td>1.46</td><td>1.78</td><td>1.80</td><td>1.66</td><td>1.42</td><td>1.48</td><td>1.50</td><td>1.35</td><td>1.38</td><td>1.42</td></tr><tr><td>1M</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>32</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td></td><td></td><td>0.85</td><td>0.86</td><td>1.06</td><td>1.23</td><td>1.22</td><td>1.48</td><td>1.05</td><td>1.04</td><td>1.38</td><td>1.04</td><td>1.04</td><td>1.04</td></tr><tr><td rowspan="4"></td><td>10K</td><td>128</td><td>96</td><td>64</td><td>128</td><td>96</td><td>128</td><td>96</td><td>96</td><td>256</td><td>128</td><td>256</td><td>128</td></tr><tr><td>100K</td><td>1.27</td><td>1.15</td><td>1.15</td><td>1.16</td><td>1.16</td><td>1.16</td><td>1.34</td><td>1.23</td><td>1.24</td><td>1.21</td><td>1.22</td><td>1.15</td></tr><tr><td>1M</td><td>64</td><td>64</td><td>64</td><td>128</td><td>64</td><td>64</td><td>128</td><td>128</td><td>64</td><td>64</td><td>128</td><td>64</td></tr><tr><td></td><td></td><td>1.18</td><td>1.09</td><td>1.13</td><td>1.12</td><td>1.10</td><td>1.11</td><td>1.25</td><td>1.21</td><td>1.28</td><td>1.10</td><td>1.12</td><td>1.10</td></tr><tr><td rowspan="4"></td><td>10K</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td><td>512</td></tr><tr><td>100K</td><td>1.21</td><td>1.27</td><td>1.15</td><td>1.15</td><td>1.16</td><td>1.16</td><td>1.34</td><td>1.23</td><td>1.24</td><td>1.21</td><td>1.22</td><td>1.15</td></tr><tr><td>1M</td><td>64</td><td>64</td><td>64</td><td>128</td><td>64</td><td>64</td><td>128</td><td>128</td><td>64</td><td>64</td><td>128</td><td>64</td></tr><tr><td></td><td></td><td>1.18</td><td>1.09</td><td>1.13</td><td>1.12</td><td>1.10</td><td>1.11</td><td>1.25</td><td>1.21</td><td>1.28</td><td>1.10</td><td>1.12</td><td>1.10</td></tr></tbody></table>

A final experiment compares the performance of our software with that of HOOMD Anderson et al (2020) on a Brownian dynamics simulation. This experiment was carried out on an A40 GPU with a Lennard-Jones interaction force, varying the total
<figure><img src="image_10.png" /><figcaption>Fig. 8 Computation times (in seconds) obtained by our software and HOOMD for a LJ simulation of 1M iterations with 10K, 100K & 1M particles and densities 5%, 10% & 20% on an A40 GPU.</figcaption></figure>

number of particles and the density over a simulation of 1 million iterations, in a
similar way to our previous tests. The results can be seen in figure 8 using a loga-
rithmic scale, which shows that our approach outperforms HOOMD by a factor of 1.2
on average. These findings demonstrate that our technique is a viable alternative at
least on an A40 GPU, even if it is challenging to compare our results with those of
more generic software like HOOMD, which can be used for many kinds of molecular
dynamics simulations and runs on both CPU or GPU architectures.

5 Implementation details

This section describes various aspects of how our code works. All files are available as
an archive supplied with the paper.

5.1 Simulation setup

Our simulation code relies on several libraries:

829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
• The CURAND library<fnref n="1" /> provides functions for generating random numbers on the GPU, which is essential for simulating Brownian motion and other stochastic processes within the particles.

• Thrust<fnref n="2" /> offers functionalities for working with device pointers, which are necessary to manage data transfer between the CPU and GPU memory during the simulation.

• The "JSON for Modern C++" library<fnref n="3" /> is used to read the configuration file and parse its contents.

• The code utilizes standard C++ libraries for various tasks such as input/output operations, data formatting, string manipulation, and random number generation on the CPU (potentially for initialization).

• Additionally, the code includes custom-made classes residing in separate header files to handle kernel optimization for efficient GPU execution (see section 3.3) and calculations related to specific interaction potentials (see the following section).

The simulation starts with a JSON file to configure various aspects of the run, an example is given at the end of the paper. The "PARTICLES" section of the JSON file specifies information for each particle type, the code iterates through this section and creates particle objects. Each object stores properties like particle type, radius, number of particles, and interaction potential parameters, allowing to define several sets of particle, each with its own parameters. The distribution for particle radii is handled by a separate class based on the "distribution" property.

The "SOLVANT" section provides properties of the surrounding solvent, including dielectric constant and viscosity. These properties are used in calculations related to particle interactions or Brownian motion. The "SIMULATION" section defines various parameters governing the simulation run. These include simulation temperature, timestep, total simulation time, output directory and saving frequency. The code uses these parameters to initialize the simulation.

<footnote n="1">¹https://docs.nvidia.com/cuda/curand/</footnote>

<footnote n="2">²https://nvidia.github.io/cccl/thrust/</footnote>

<footnote n="3">³https://github.com/nlohmann/json</footnote>
Finally, the "ACCELERATION" section gives access to the "Rc_multiplier" 921
parameter described in section 3.2 to control the update frequency of neighbour lists 922
during the simulation, and the "maxListSize" parameter gives the maximal size for 923
these lists. 924

## 5.2 Definition of interaction forces

As explained earlier, the code corresponding to Lennard-Jones, Yukawa and DLVO 931
interaction forces is directly available in a separate directory. The "FORCE" section in 932
the JSON file allows to choose between these interaction forces and specify parameters 933
such as the potential depth or the cutoff radius. 934

Based on the interaction forces already defined, a new type of force can be added 938
by applying the following steps: 939

* A C++ class must be defined with a structure describing the parameters used by the 941
  potential. The exact meaning of each parameter depends on the specific potential 942
type, but generically they might represent the cutoff radius, the potential strength, 943
etc. 944

* The "init" function handles reading parameters values from the JSON file. 948

* The "compute" function calculates the force acting between two particles based on 950
  their positions (usually denoted "r0" and "r1") and other specific parameters. 951

* The section of code responsible for calculating the total force acting on each particle 953
  must be modified to include a call to the new force calculation function. 954

* If the force depends on particle properties beyond mass and position (e.g., charge), it 957
  might also be necessary to modify particle data structures to store these additional 958
properties. 959

If the new force shares certain characteristics with existing forces, an experienced 962
programmer may also consider creating a new class using inheritance mechanisms to 963
implement the new force's behavior. 964
## 5.3 Restart procedure

Like most MD simulation software, our code offers flexibility in saving simulation results for visualization and analysis purposes. It allows to export particle data at regular intervals or logarithmically spaced intervals if capturing data more frequently is necessary at the beginning of the simulation. The exported data is stored in a modified XYZ file format which includes the spatial coordinates of each particle in a simple text-based structure with a user-controlled precision. Several software applications, including Ovito Stukowski (2010), can read such a file to visualize and analyze the trajectory of simulated particles.

A binary file storing double-precision particle positions is also created after each save in XYZ format. Only one copy of this binary file is kept, i.e. only the positions calculated during the last computation iteration are accessible in double precision. This feature was made necessary by the fact that the GPUs available in a computing cluster are sometimes subject to a time limit, which can abruptly interrupt a simulation in progress. Our system thus guarantees that the simulation can be restarted on any machine from the last binary save without introducing position-related calculation errors, which would not be the case with an XYZ file. This also makes it possible to restart the simulation from positions calculated with other software.

## 5.4 Potential improvements

The use of Thrust for sorting introduces a slight delay due to memory allocation, as it doesn't perform in-place sorting by default. Although custom allocators for Thrust sort could indeed provide performance improvements, they were not implemented in our code due to the need for extensive testing. Initial experiments indicated a modest gain, with sorting times reduced by 2-3 ms for a 32M array on a 3090 GPU (from 140 ms with the default Thrust implementation to 137 ms with custom allocation). While this
suggests potential benefits, we prioritize stability and more impactful optimizations for now.

The fact that our code is limited to double precision raises interesting questions about the trade-off between numerical accuracy and computing efficiency. Switching to single precision, while potentially improving computational efficiency on GPUs, would introduce unacceptable divergence in the final results due to accumulation errors and compromise the reliability of the simulations. An alternative approach, as used in other software like HOOMD, would be to adopt reduced units, which can mitigate precision issues while maintaining efficiency. Exploring the use of reduced units could be a promising solution for the future, allowing for both accurate results and improved computational performance.

# 6 Conclusion

In this paper, we presented a high-performance BD simulation code designed for efficient execution on modern Nvidia GPUs. The code leverages key features like double-precision calculations, a bucket sort neighbor list management strategy, and a CUDA occupancy-aware autotuning system to achieve optimal performance across various GPU architectures. The software offers flexibility in simulating particles interacting through Lennard-Jones, Yukawa, and DLVO potentials, with the capability to incorporate additional interaction forces through a well-defined extension mechanism.

Our benchmarks demonstrate significant performance improvements compared to existing tools when dealing with large-scale BD simulations. The code offers user-friendly configuration options via a JSON file and functionalities like restarting simulations from both binary files and XYZ data, ensuring smooth workflow integration.
1059 Further research and development efforts will focus on expanding the supported
1060 interaction potential library, including additional interaction forces commonly encoun-
1061 tered in material science simulations. We would also like to refine the autotuning
1062 system to adapt to a wider range of hardware architectures and complex interaction
1063 potentials, and conduct new benchmarks on high-performance computing clusters with
1064 multiple GPUs to assess the software's scalability for large-scale simulations.
1065
1066
1067
1068

1069 We believe this software provides a valuable tool for researchers in material science
1070 and related fields, enabling efficient and accurate BD simulations. The open-source
1071 nature of the code facilitates further development and customization by the scientific
1072
1073
1074 community.
1075
1076

1077
1078

## Acknowledgements

1079 This research was funded, in whole or in part, by the French National Research Agency
1080 (ANR), project ANR-20-CE46-0004. Authors thank also Region Nouvelle-Aquitaine
1081
1082 for the cofunding of the SOMA-DNS project.
1083
1084
1085

1086
1087

## References

1088 Abraham MJ, Murtola T, Schulz R, et al (2015) GROMACS: High performance molec-
1089 ular simulations through multi-level parallelism from laptops to supercomputers.
1090
1091 SoftwareX 1-2:19-25
1092
1093
1094 Adjoua O, Lagardère L, Jolly LH, et al (2021a) Tinker-HP: Accelerating molecular
1095 dynamics simulations of large complex systems with advanced point dipole polariz-
1096
1097 able force fields using GPUs and multi-GPU systems. Journal of Chemical Theory
1098
1099 and Computation 17(4):2034-2053
1100
1101
1102
1103
1104
<table>
  <tbody>
    <tr>
      <td>Adjoua O, Lagardère L, Jolly LH, et al (2021b) Tinker-HP: Accelerating molecular dynamics simulations of large complex systems with advanced point dipole polarizable force fields using GPUs and multi-GPU systems. Journal of Chemical Theory and Computation 17(4):2034–2053. PMID: 33755446</td>
      <td>1105<br>1106<br>1107<br>1108<br>1109<br>1110<br>1111</td>
    </tr>
    <tr>
      <td>Anderson JA, Glaser J, Glotzer SC (2020) Hoomd-blue: A python package for high-performance molecular dynamics and hard particle monte carlo simulations. Computational Materials Science 173:109,363</td>
      <td>1112<br>1113<br>1114<br>1115<br>1116<br>1117</td>
    </tr>
    <tr>
      <td>Bender, Jan (2024) https://github.com/InteractiveComputerGraphics/cuNSearch</td>
      <td>1118<br>1119<br>1120</td>
    </tr>
    <tr>
      <td>Cerbelaud M, Tran CT, Ferrando R, et al (2017) Interdiffusion and crystallization of oppositely charged colloids. Phys Chem Chem Phys 19:31,094–31,102</td>
      <td>1121<br>1122<br>1123<br>1124</td>
    </tr>
    <tr>
      <td>Chisholm R, Maddock S, Richmond P (2020) Improved GPU near neighbours performance for multi-agent simulations. Journal of Parallel and Distributed Computing 137:53–64</td>
      <td>1125<br>1126<br>1127<br>1128<br>1129<br>1130</td>
    </tr>
    <tr>
      <td>Colberg PH, Höfling F (2011) Highly accelerated simulations of glassy dynamics using GPUs: Caveats on limited floating-point precision. Computer Physics Communications 182(5):1120–1129</td>
      <td>1131<br>1132<br>1133<br>1134<br>1135<br>1136</td>
    </tr>
    <tr>
      <td>Dice K, Gray C, Walker G, et al (2018) CUDA-accelerated simulation of Brownian dynamics. In: Proceedings of the Practice and Experience on Advanced Research Computing. Association for Computing Machinery, PEARC ’18</td>
      <td>1137<br>1138<br>1139<br>1140<br>1141</td>
    </tr>
    <tr>
      <td>Fernández-Fernández JA, Westhofen L, Löschner F, et al (2022) Fast octree neighborhood search for SPH simulations. ACM Trans Graph 41(6)</td>
      <td>1142<br>1143<br>1144<br>1145<br>1146</td>
    </tr>
    <tr>
      <td>Gross J, Köster M, Krüger A (2019) Fast and Efficient Nearest Neighbor Search for Particle Simulations. In: Computer Graphics and Visual Computing (CGVC)</td>
      <td>1147<br>1148<br>1149<br>1150</td>
    </tr>
  </tbody>
</table>
1151 Hashemi A, Peláez RP, Natesh S, et al (2023) Computing hydrodynamic interactions
1152 in confined doubly periodic geometries in linear time. The Journal of Chemical
1153 Physics 158(15):154,101
1154
1155
1156
1157 Hijma P, Heldens S, Sclocco A, et al (2023) Optimization techniques for GPU
1158 programming. ACM Comput Surv 55(11)
1159
1160
1161 Hoetzlein R (2014) Fast and Efficient Nearest Neighbor Search for Particle Simula-
1162 tions. In: GPU Technology Conference (GTC)
1163
1164
1165 Howard MP, Anderson JA, Nikoubashman A, et al (2016) Efficient neighbor list calcu-
1166 lation for molecular simulation of colloidal systems using graphics processing units.
1167
1168 Computer Physics Communications 203:45-52
1169
1170
1171 Howard MP, Statt A, Madutsa F, et al (2019) Quantized bounding volume hierar-
1172 chies for neighbor search in molecular simulations on graphics processing units.
1173
1174 Computational Materials Science 164:139-146
1175
1176
1177 Israelachvili JN (2011) Intermolecular and Surface Forces. Academic Press
1178
1179 Kim H, Fábián B, Hummer G (2023) Neighbor list artifacts in molecular dynamics
1180 simulations. Journal of Chemical Theory and Computation 19(23):8919-8929
1181
1182
1183 Lipscomb TJ, Zou A, Cho SS (2012) Parallel Verlet neighbor list algorithm for
1184 GPU-optimized MD simulations. In: Proceedings of the ACM Conference on Bioin-
1185 formatics, Computational Biology and Biomedicine. Association for Computing
1186
1187
1188 Machinery, New York, NY, USA, BCB '12, p 321-328
1189
1190
1191 Ohno K, Nitta T, Nakai H (2017) SPH-based fluid simulation on GPU using Ver-
1192 let list and subdivided cell-linked list. In: 2017 Fifth International Symposium on
1193
1194 Computing and Networking (CANDAR), pp 132-138
1195
1196
Stukowski A (2010) Visualization and analysis of atomistic simulation data with OVITO - the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering 18(1)

Thompson AP, Aktulga HM, Berger R, et al (2022) LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications 271:108,171

Tran CT, Crespin B, Cerbelaud M, et al (2018) Colloidal suspension by SRD-MD simulation on GPU. Computer Physics Communications 232:35-45

Winkler D, Rezavand M, Rauch W (2018) Neighbour lists for smoothed particle hydrodynamics on GPUs. Computer Physics Communications 225:140-148

Zhu Y (2022) RTNN: Accelerating neighbor search using hardware ray tracing. In: Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, ACM

## Appendix A Appendix: example of JSON configuration file

```json
{
  "PARTICLES": [
    {
      "type": 0,
      "radius": 3e-07,
      "maxRadius": 0.0,
      "minRadius": 0.0,
      "radius_std_dev": 0.0,
      "distribution": "lognormal",
    }
  ]
}
```
1243 "number": 50000,
1244 "volumicMass": 2200.0,
1245
1246 "psi":1
1247 },
1248
1249 {
1250
1251 "type": 1,
1252 "radius": 3e-07,
1253
1254 "maxRadius": 0.0,
1255
1256 "minRadius": 0.0,
1257
1258 "radius_std_dev": 0.0,
1259 "distribution": "lognormal",
1260
1261 "number": 50000,
1262
1263 "volumicMass": 2200.0,
1264
1265 "psi":-1
1266 }
1267
1268 ],
1269 "FORCE": [
1270
1271 {
1272
1273 "name": "Yukawa",
1274
1275 "k":1e8,
1276 "U":10,
1277
1278 "rcrep":1,
1279 "coeffdir":7.93e-13,
1280
1281 "Rc": 3.2,
1282
1283 "first": -1,
1284
1285 "second": -1
1286 }
1287
1288
1289
],
1290
"SOLVANT": {
1291
"dielectricConstant": 81.0,
1292
"eta": 0.001
1293
1294
},
1295
"ACCELERATION":
1296
{
1297
"Rc_multiplicator": 1.05,
1298
"maxListSize": 128
1300
1301
1302
1303
},
1304
"SIMULATION": {
1305
"temperature": 293.0,
1306
"overlapDistance": 1.1,
1307
"deltaTime": 3.685e-07,
1308
"eta": 0.001,
1309
"totalTime": 0.3685,
1310
"outputDir": "out2/",
1311
"outputScale": 3333333.3333333,
1312
"savePeriod": -1,
1313
"savePrecision": 4,
1314
"restartSaveInterval": 100000,
1315
"randomSeedPosition": 1234,
1316
"density": 0.20,
1317
"restartFile":"
1318
}
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
}
1329
1330
}
1331
1332
1333
1334