GPU parallel simulation algorithm of Brownian particles with excluded volume using Delaunay triangulations

TL;DR

This paper introduces a GPU-based parallel algorithm for simulating Brownian particles with excluded volume, utilizing Delaunay triangulations to efficiently handle overlaps and achieve significant performance improvements.

Contribution

The paper presents a novel GPU parallel simulation method using Delaunay triangulations to manage particle overlaps in Brownian motion simulations, enhancing efficiency.

Findings

Achieved up to 100x speedup over sequential methods.

Validated the algorithm with active colloidal particle models.

Maintained triangulation validity during particle displacements.

Abstract

A novel parallel simulation algorithm on the GPU, implemented in CUDA and C++, is presented for the simulation of Brownian particles that display excluded volume repulsion and interact with long and short range forces. When an explicit Euler-Maruyama integration step is performed to take into account the pairwise forces and Brownian motion, particle overlaps can appear. The excluded volume property brings up the need for correcting these overlaps as they happen, since predicting them is not feasible due to the random displacement of Brownian particles. The proposed solution handles, at each time step, a Delaunay triangulation of the particle positions because it allows us to efficiently solve overlaps between particles by checking just their neighborhood. The algorithm starts by generating a Delaunay triangulation of the particle initial positions on CPU, but after that the triangulation is always kept on GPU memory. We used a parallel edge-flip implementation to keep the triangulation updated during each time step, checking previously that the triangulation was not rendered invalid due to the particle displacements. The algorithm is validated with two models of active colloidal particles. Upon testing the parallel implementation of a long range forces simulation, the results show a performance improvement of up to two orders of magnitude when compared to the previously existing sequential solution.

Equations

$${\vec{r}}_i(t+\Delta t)={\vec{r}}_i(t)+{\vec{F}}_i(t)\Delta t+\vec{\xi }\sqrt{D\Delta t},$$

$${\vec{F}}_i=\sum _{k\ne i}{\mu }_i{\alpha }_k\vec{f}({\vec{r}}_i-{\vec{r}}_k),$$

$${\vec{F}}_{ij}({\vec{r}}_{ij})=\{\begin{array}{ll}{\mu }_i{\alpha }_k\vec{f}({\vec{r}}_{ij}),& \text{if }{r}_{ij}\le r_{\text{cutoff}}\\ 0,& \text{otherwise}\end{array}$$

$${\vec{r}}_1\prime$$

$${\vec{e}}_0$$

$${\vec{e}}_1$$

$${\vec{e}}_2$$

$$\text{flip}({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3,{\mathbf{v}}_4)=\{\begin{array}{ll}(s_2\le 0)\wedge (t_2\le 0)\wedge (s_2+t_2\ge d_2),& \text{if }{d}_2<0\\ (s_2\ge 0)\wedge (t_2\ge 0)\wedge (s_2+t_2\le d_2),& \text{otherwise}\end{array}$$

$$L(N)$$

$${\vec{r}}_i(t+\Delta t)={\vec{r}}_i(t)+V_0(\cos {\theta }_i\hat{x}+\sin {\theta }_i\hat{y})\Delta t.$$

$${\theta }_i(t+\Delta t)={\theta }_i(t)+\sqrt{2D\Delta t}{n}_i,$$

Full text

††journal: Computer Physics Comunnications

1 Introduction

A colloidal suspension is a mixture of microscopical insoluble particles dispersed throughout a continuous fluid, where particle sizes range from 1 nm to 10 $\mu$m. Colloidal suspensions appear in several natural and artificial substances as the milk, mud, inks, cosmetics or latex paint, for example. Also, they are used in many intermediate industrial processes The interactions between colloidal particles of various kinds [1] have effects on the physical and chemical properties of the mixture such as its viscosity or light dispersion. To study these and other properties it is necessary to simulate particle systems of growing numbers ($N\geq 10^{4}$). Also, colloids are being used as models for active systems, to describe the motion of self-propelled microorganisms [2, 3, 4].

Colloids can be modelled as hard bodies subject to Brownian diffusive motion. Colloidal particles can typically interact through the fluid in what is called hydrodynamic interactions, via electrostatic forces for charged colloids, which can be screened in an electrolyte, or with van der Waals forces [1]. In out of equilibrium conditions, phoretic forces also appear [5]. Except for the hydrodynamic forces, these interactions can be modelled with good approximation as pairwise additive forces, which in out of equilibrium conditions can eventually break the action-reaction symmetry.

The simulation of colloidal dispersions, can be divided on two main problems executed in sequence: updating the positions of the particles due to the interparticle interactions, according to some integration rule and ensuring that the bodies do not overlap because of their movement, in order to respect the excluded volume interaction. These problems are specific instances of the n-body simulation and collision detection respectively [6]. In some contexts, the simulation of colloidal particles is referred as Brownian dynamics.

There are two main methods for solving overlaps between particles: correcting all of them at once after they happen or use an event-driven approach, integrating the system until the collision instant, process the involved particles and repeat until the system reaches the target time step. The last method, is particularly useful when inertia is important and collisions result in rebounces as in granular materials [7, 8]. It requires knowing the positions of the involved bodies at the time of collision, which becomes difficult when random Brownian motion is present. For the simulation of colloidal particles, which lack of inertia and excluded volume acts like a boundary condition rather than producing collisions, the first method is more suited.

This work focuses on designing and implementing a novel parallel simulation algorithm for 2D colloidal particle interacting with short and long range pairwise forces, with periodical boundaries, excluded volume and Brownian motion. The algorithm implementation takes advantage of the data-parallel computing capabilities of the GPU architecture, which have proven to be effective at accelerating the simulation process of several computational physics problems [9, 10, 11]. The interactions forces are allowed to be non-reciprocal as in the case of active particles [12, 13]. The main contribution of this work consists of a new and efficient method of resolving particle overlaps by using Delaunay triangulations, which are maintained periodically and fully on the graphics card. The starting positions and triangulation are initialized on the host while all the simulation code is executed on the device. The random values are also generated on the graphics card, both on the initialization and simulation phases. The algorithm uses a GPU edge-flip implementation to keep the triangulation fulfilling the Delaunay condition during each time step and to correct inverted triangles in case they are generated due to the particle displacements. For the short range force simulation, we developed a parallel algorithm that builds and uses Verlet lists in order to handle the particle neighborhood in parallel. The algorithm is validated with two models of active colloidal particles. Upon testing the parallel implementation of a long range forces simulation, the results show a performance improvement of up to two orders of magnitude when compared to the previously existing sequential solution. The algorithm for the short range force presents a similar performance improvement regarding the parallel long range implementation.

The paper is organized as follows: Section 2 describes the specific conditions and properties that the simulated systems must operate under. Section 3 lists previous related work used to solve similar problems. Section 4 details the designed solution with its subcomponents, data structures, and optimizations. The implementation of the algorithm is described in Section 5. Sections 6 and 7 cover the tests, benchmarks, validation and used methodology, presenting the running time and performance results when compared to the other implemented solutions. Finally, section 8 rounds up the obtained results.

2 Desctiption of the model

This section contains the description of the involved concepts and properties of this problem that may differentiate it from other body simulation problems, such as the excluded volume and stochastic component of particle movement.

2.1 Preliminaries

Let $P=\{p_{1},p_{2},...,p_{N}\}$ be a set of $N$ bodies on a d-dimensional space. The n-body simulation is the computation of the interactions over each body in $P$, where $F_{i}$ corresponds to the interaction over $p_{i}$ by effect of $P_{i}=P\setminus\{p_{i}\}$, the set of all other bodies in the system. The interaction forces $F$ typically depend on the distance $r_{ij}$ between two bodies as $F\sim r_{ij}^{-q}$. If $q>d$ the force is said to be short ranged, while if $q\leq d$, it is considered a long range force. When the forces are long ranged, the set $P_{i}$ cannot be reduced in n-body simulations and an exact evaluation of the forces has a cost $O(N^{2})$. Approximate solutions for long range interactions as the Barnes-Hut algorithm reduce the cost to $O(N\log N)$ [14]. But for short range forces, the interactions can be truncated and, therefore, $P_{i}$ can be reduced to the neighborhood of particles close to $p_{i}$. In this case, the evaluation of the forces costs $O(N*N_{NL})$ on average, where $N_{NL}$ is the average number of neighbors of a body [15]. Since the construction of the list is $O(N^{2})$ for evaluating all pairwise distances between bodies, it is possible to partition the simulation domain in cells so that close bodies get binned together in the same cells. Assignment of bodies to their respective cells takes $O(N)$ time [16, 17].

The simulation domain is a two-dimensional $L\times L$, periodical box across the $X$ and $Y$ axes, meaning that the particles wrap around the box as they move across its boundaries. For the distance calculations between particles, including force calculations, we follow the minimum-image convention, in which a particle interacts with another via its real position or its image depending on which is the shortest.

2.2 Particle interaction without excluded volume

Microscopic particles move in an overdamped regime, with no inertia. When subject to a force $\vec{F}$, the equation of motion is simplified to $d\vec{r}/dt=\gamma\vec{F}$, where $\gamma$ is the mobility. Absorbing the mobility coefficient into the force, which will then have velocity units, in a time step $\Delta t$, the integration rule for updating the position of a particle over time is performed using the Euler-Maruyama method:

[TABLE]

where $\vec{F}_{i}$ is the deterministic velocity obtained from the interactions between the particle $i$ and $P_{i}$, $D$ is the diffusion coefficient, and $\vec{\xi}$ is a random vector, where the components follow a normal distribution of zero mean and unit variance, and corresponds to a noise added that takes into account the diffusive Brownian motion [18].

The force model we use for the simulations describes the interaction of self-diffusiophoretic active particles [12]. In this model, particles can be of different type, characterized by two charges, $\alpha$ and $\mu$; the former is responsible of creating the concentration field, while the second describes the response of a particle to the field, leading to the following interaction law:

[TABLE]

where $\vec{f}(\vec{r})={\vec{r}}/{r^{3}}$ for the studied long range force, while $\vec{f}(\vec{r})={\vec{r}}/{r^{7}}$ for the short range interaction. Note that if $\alpha_{i}\neq\mu_{i}$, the action-reaction symmetry is broken and self-motion is possible. Charged colloidal particles are included in this model if $\alpha_{i}=\mu_{i}=q_{i}$, equal to the electric charge of the particles.

Since the short range force decays much faster with distance compared to the long range force, its calculation considers a cutoff radius from which the value of the force is considered zero, as shown on Figure 1. The short range force is then computed as:

[TABLE]

We used $r_{\text{cutoff}}=2.5\sigma$ for the simulated short range force in our experiments, where $\sigma$ is the particle diameter.

2.3 Excluded volume

The simulated particles are represented as hard disks with a uniform diameter $\sigma$. Although here we consider only monodisperse colloids, it is direct to extend the method to polydisperse systems where radii dot not differ too much. Since the integration rule (1) ignores the excluded volume condition, it can happen that the updated positions produce overlaps between two or more particles, resulting in a physical impossibility. In order to ensure this property, at the end of each time step, the members of all overlapping pairs $(p_{i},p_{j})$ are moved apart from each other in a way that corrects the overlaps:

[TABLE]

where $\vec{r}_{1}$ and $\vec{r}_{2}$ are the original positions and $\vec{r}_{1}\prime,\vec{r}_{2}\prime$ the updated positions. If $\delta^{*}=(\sigma-|\vec{r_{12}}|)/2$, the particles would move in opposite directions from each other along $\hat{n}=(\vec{r}_{2}-\vec{r}_{1})/|\vec{r}_{12}|$, leaving the particles in tangential contact. If $\delta^{*}=\sigma-|\vec{r}_{12}|$, the movement is proportional to the magnitude of the previously existing overlap, simulating a bounce effect resulting from the collision at some instant $t^{*}\leq t+\Delta t$. This last value is the one used for processing the overlaps in the simulation and guarantees that no accumulation is produced at the contact distance.

2.4 Stochastic displacements

The particle displacements on (1) have a random noise component $\vec{\xi}$, modeled as a random variable with standard normal (or Gaussian) distribution with zero mean and standard deviation $\hat{\sigma}=\sqrt{D\Delta t}$. Reducing $\Delta t$, the deterministic and stochastic displacements in each time step are reduced. However, for a Gaussian distribution, it is always possible that large values are generated (at the tail of the distribution), leading to excessively large displacements (see Figure 2). To avoid these problems, the simulation ignores values larger than 3 standard deviations. We considered two methods in order to achieve this, as shown on Figure 3:

1. (a)

Reroll the values outside the range $[-3\hat{\sigma},3\hat{\sigma}]$. 2. (b)

Truncate the values to the range $[-3\hat{\sigma},3\hat{\sigma}]$.

Both methods produce probability distributions different from each other and from the original; while the first alternative raises the probability of all values in range, the second one raises the probability at the edges. These modifications do not generate a noticeable statistic distortion, since the considered range includes $99.7\%$ of the possible values. In our simulations, we opted for the second method, which turns out to be faster and better suited for parallel execution, since it needs to generate a single random number instead of a variable quantity of random values in the first method.

3 Related work

For short-range forces calculation, the standard technique is the use of Verlet lists [19, 16, 17]. The authors in [20, 21] parallelize the list construction by having a $O(N^{2})$ list of all possible pairs of bodies. A predicate checking closeness between the pair members is evaluated over all elements of the list, which can then be used for a key-value sort to group all the neighboring pairs consecutively in the array. A parallel scan operation allows to get the number of elements that must be copied to the neighbor list. The authors then combine this algorithm with fixed cell partitioning in order to replace distance calculations with less-expensive cell neighborhood checks.

For the parallel n-body simulation, with full calculation of the $O(N^{2})$ forces, Nyland et al. [22] developed a grid-style tiling algorithm, reading the particles from the global space and storing them on GPU shared memory, increasing performance as multiple threads read from that space at a lower latency. Partitioning the load/store process on groups of $p$ particles allows fitting an arbitrary input size on the hardware-limited shared memory size. Burtscher and Pingali [23] parallelized the Barnes-Hut simulation [14], which computes an approximation for the force, representing the cell hierarchy k$d$-tree as multiple arrays for each node field. It uses atomic lock operations to build the tree in parallel, throttling the threads that failed to get the lock so they do not waste bandwidth with unsuccessful lock requests. The tree is then filled with the center of mass data, starting from lower nodes in the tree according to the order of allocation for the scan. Bedorf et al. [11] uses a Z-order curve to sort the particles spatially. Each thread is assigned to a particle, applying a mask value to it to determine the octree cell the particle should be assigned. The linking of the tree is made by assigning a thread to each cell node and then doing a binary search over the corresponding Z-order key to find both the parent and child nodes, if appropiate.

To detect and process collisions, Hawick and Playne [24] developed a multi-GPU algorithm with a tiling scheme similar to the one used by Nyland et al. [22]. If a pair of particles overlap, the associated threads store the index of its colliding neighbor and the time at which the collision occurred. The collisions then are resolved iteratively starting from the earliest, redoing the previous process in order to find possible new collisions.

Finally, for overlap correction, Strating [25] describes a brute-force sequential algorithm that checks all pairs of bodies for possible overlaps and corrects them following equation (4). The algorithm may need to iterate an unbounded number of times at each time step because some corrections may generate new overlaps with neighboring particles.

4 Algorithm

This section describes the parallel algorithm in detail. It includes the generation of the initial data, data structures, overlapping detection and correction, and Delaunay condition updates, among others, putting emphasis in what threads are doing at each time step.

4.1 Overview

The simulation consists of two phases: (i) sequential initialization of the simulation data, followed by a host to device transfer and (ii) a parallel simulation phase. The initial positions are initialized over a triangular mesh with $N^{*}\geq N$ vertex, where each vertex represents a particle and their types are assigned randomly according to the specified concentrations. A sample of $N$ particles is selected from the mesh by Reservoir Sampling [26], resulting in the input particle set, which is homogeneous in space.

The n-body algorithm for the long range force is based on a grid-style tiling, which uses the shared memory of the multiprocessor assigned to each thread block to store the particles in groups. In this algorithm each thread is mapped to exactly one particle in the system, and since it is possible to lack action-reaction symmetry on the force, no redundant computation is done unlike the cases where $\vec{f}_{ij}=-\vec{f}_{ji}$.

Once the forces are calculated, the particles are advanced one time step using eqn (1). As a result, particle overlaps can appear. When $\Delta t$ is small enough, for hard disks of similar or equal radii, only neighbor particles can overlap. Then, to detect and correct overlaps, instead of a brute-force algorithm that would check all $O(N^{2})$ pairs, only neighbors are checked. The Delaunay triangulation [27] is particularly well suited to detect neighbors for monodisperse or slightly polydisperse disks. In dense systems, the overlap corrections can be highly non-local, as the correction of on pair can generate a sequence of other overlaps that need correction. It is therefore not clear a priori the computational cost of this stage, which is the reason why we consider both short and long range interaction forces.

The Delaunay triangulation can be built constructively or from an existing triangulation. Lawson’s algorithm [28] accepts a triangulation as input and transforms it into a Delaunay triangulation via a finite sequence of edge-flip operations [27]. Based on the Lawson algorithm, Navarro et al. [29] developed a parallel implementation for generating quasi-Delaunay triangulations, so it is possible to keep the triangulation updated without need of host-device memory transfers. These are quasi-Delaunay because exact predicates are too expensive on the GPU; nevertheless, this approximate construction is sufficient for our problem. However, since the input for the Navarro et al. algorithm must be a valid triangulation, we must first correct potential triangle inversions with invalid edge intersections, which can result from the particle displacements. With this strategy, the Delaunay triangulation of the particle positions is built only once from scratch on the host [30, 31], which is transferred to the device. Thereafter, the triangulation is maintained updated after each time step on the device.

4.2 Data structures

We store the simulation data as a Structure of Arrays (SoA) on global device memory, using total $O(N)$ space. The particle data consists of their position $(x_{i},y_{i})$ and their charges $(\alpha_{i},\mu_{i})$, stored as floating point vector types 111float2 or double2. in order to increase bandwidth utilization [32]. We use an additional buffer array for positions so that writes are not done at the same adresses for reads, avoiding a synchronization step. We store the simulation parameters that remain unchanged during a same instance on a constant device memory structure [32], such as $N,D,\Delta t,\sigma$ and derived constants $\sigma^{2}$ and $\sqrt{D\Delta t}$. Additionally, we store the triangulation data using the same scheme as [29].

4.3 Inverted triangle detection

There are two possible reasons for a triangle inversion: an edge gets inverted because the distance vector between its terminal vertex changes orientation or because one of its vertex crosses the edge opposite to it. The first problem means that the particles went through each other, which is a physical impossibility and must not be allowed, and can only takes place for large $\Delta t$. The criterion $\vec{r}_{ij}^{0}\cdot\vec{r}_{ij}^{1}<0$ is used to determine if the above situation happened on the current time step, where $\vec{r}_{ij}=\vec{r}_{j}-\vec{r}_{i}$ is the distance vector between the compared particles, and $\vec{r}_{ij}^{0}$ and $\vec{r}_{ij}^{1}$ are evaluated with current and previous positions, respectively. The buffer array allows to compare the distances before and after integration, and the edges of the triangulation show the pairs that need checking. Once an invalid movement is detected, the last positions are discarded and integration is repeated with a lower $\Delta t$ value than the currently used.

Finally, the inverted triangle detection becomes equal to checking if a vertex crossed (towards a neighboring triangle) any of the edges that enclose it, which is equivalent to point-in-triangle detection. Using the barycentric coordinates $d$, $s$, and $t$ on triangles (see Figure 4),

[TABLE]

the criterion used to detect the edges that must be flipped is

[TABLE]

Figure 4 shows the vectors used on the predicate that checks if the original edge must be flipped. The predicate becomes true when applied to edge $(v_{1},v_{2})$ in (4(b)), because point $v_{4}$ lies inside the triangle $(v_{1},v_{2},v_{3})$. This means that $(v_{1},v_{2})$ must be replaced with $(v3_{1},v_{4})$, as in a common edge-flip operation. In this case, $v_{3}$ and $v_{4}$ are stored as opposite vertices to edge $(v_{1},v_{2})$ in the triangulation data structure in GPU device memory.

It is worth noting on Figure 5 that the movement of point $b$ across edge $(c,d)$ creates an intersection between it and edge $(b,e)$. The edge flip between $(a,b)$ and $(c,d)$ removes the inverted triangle, restoring the local triangulation. The triangulation may still not satisfy the Delaunay property, so additional edge flips may be needed on further steps. For this stage, we use the Navarro et al. algorithm [29].

Figure 6 shows a situation where particle $f$ moves across edges $(b,e)$ and $(c,e)$, needing two consecutive flips in order to restore the local triangulation. While function 8 cannot properly evaluate this case or similar movements across further distances, this kind of inversion can be detected allowing to revert the step. Anyway, small time steps guarantee that this situation is extremely unlikely to happen.

4.4 Overlap correction

The overlap correction uses the topological information contained in the edges of the Delaunay triangulation, which allows for each particle fast access to the neighborhood of particles that may be overlapping with it. The algorithm maps threads to edges in such a way that each thread handles one edge of the triangulation. A thread gets the positions of the particles that form the edge, checking if there exists an overlap between them ($r_{ij}<\sigma$). If the check is positive, the algorithm computes the displacements of the involved particles according to (3) Since the same particle can be part of many edges at once, the algorithm sums atomically the displacements in a global array, in order to avoid concurrency hazards.

Once the algorithm computes the total displacements, it maps each thread with a particle in the same way as described previously for edges. Each thread then updates the position of its particle, applying the periodic boundary conditions when necessary. It is possible that the updated positions may still have some of the previous overlaps or even have some newly generated ones. In this case, the algorithm repeats the previous process until no overlaps are present.

When adding the partial displacements on a particle, it may happen, if these point in the same direction, that the resulting total displacement is excessively large (see Figure 7). These displacements, larger to what is needed to solve the overlap, can generate new overlaps. Eventually, the correction of the new overlaps can result in an instability, where the displacements increase with alternating sign and the iterative procedure does not converge. A solution is truncating the displacement with the heuristic value $\sigma/4$ (half the particle radius), preventing the emergence of the instability.

Finally, it is worth mentioning that the parallel correction algorithm presented here does not correspond to a parallelization of the sequential algorithm of Strating [25], which displaces particles sequentially, while in our case the displacements are added and performed in parallel. Hence, due to the chaotic dynamics of the system, the small differences in these algorithms will produce different outputs for finite $\Delta t$.

4.5 Long range forces

An improvement to the long range force calculation consists on using the intrinsic warp shuffle instruction, which allows a thread access to the registers of other threads belonging to the same warp. Each thread is assigned a lane number that identifies it from the other warp members, allowing them to read different particles from global memory. Then, each warp member takes turns in propagating the data of its corresponding particle to the other threads, who can read it via the _shfl() instruction by passing as argument the lane number of the thread currently in turn. Once the whole warp has shared the data among its members, each member reads a particle from global memory and repeats the same process until all particles have been visited. The main advantage of this optimization is a greater efficiency of memory accesses, since most of the time the threads are sharing data at registry level instead of more expensive load requests on global memory. Also, the concurrent execution of the warp members makes unnecessary the explicit synchronization of the inner warp shuffle loop.

The $.x$ and $.y$ operators reference the data of the vectorized CUDA structures for each respective variable. For example, the noise $\vec{\xi_{i}}$ has $\hat{x}$ and $\hat{y}$ components, so it is grouped as a single vector for increased memory performance [32].

5 Implementation

The parallel algorithms described in the previous section were implemented on CUDA 7.5 and C++ 11, using function templates to choose between float and double precision formats at compile time. We use the CGAL library [30] to create the 2D periodic Delaunay triangulation, which is then sent to device memory alongside the particle data before starting the simulation. The random numbers used on the parallel implementations when initializing the starting positions and generating noise during integration are created with the XORWOW pseudorandom number generation algorithm of the cuRAND library [33], using the host and device APIs respectively. Each configuration has two particle types, although the program can support a variable number of particle types for simulating. For comparison purposes, we also implemented a fully sequential long range forces algorithm with the overlap correction discussed on [25], and a parallel short range forces algorithm using Verlet lists and a discrete grid over the simulation box. The neighbor list computation during the Verlet lists construction in parallel is similar to the one described on [34], grouping together all the particles that belong in the same cells.

When the simulation finishes, the final positions are brought back to host memory and written to an output file. The visualizations on Figures 8, 13, 14 and 15 were generated reading the respective output files.

6 Performance results

Parameters:

We generated inputs for 11 different values of $N$ and 5 parameter configurations, as described on Table 1. The starting positions are generated semi-randomly as described on section 4, keeping the same seed value for the random number generator across all simulation instances.

Each configuration has two types of particles with charges $\alpha_{i},\mu_{i}$. The fraction of particles of each type is given by $\phi_{i}=N_{i}/N$, where $N_{i}$ is the number of particles of each type and $N=N_{1}+N_{2}$ is the total number of particles. The area fraction $\rho=N(\sigma/2)^{2}/L^{2}$ is a measure of the particle density. To study scaling times, we change the length of the simulation box to keep density constant when increasing $N$:

[TABLE]

The values for $N$ start at $2^{10}$, raising the exponent by 1 until $N=2^{20}$. Finally, we kept constant the values for $\sigma=1,\Delta t=0.01,D=0.01,\delta=1.0$ for all configurations and input sizes.

Figure 8 displays the particle positions after $10^{4}$ time steps for each configuration. The election of the parameters help to test the algorithms under different conditions of fluidity, density and homogeneity. For c0, there is an asymmetric attraction between particles of type 1 and 2, resulting in an homogeneous mixture, with fluid-like motion. In c1, there is a larger concentration of type 2 particles, which self-attract forming a dense cluster, segregated from type 1 particles, which self-repel forming in a gas-like state. In configuration c2, equal particles repel, while dissimilar particles attract, favouring the formation of chain-like structures, where 1 and 2 particles alternate. In c3, the situation is the opposite, where equal particles attract, while dissimilar particles repel, leading to the formation of dense segregated clusters. Finally, the interactions in c4 are analogous to those of c2, in a dilute regime, resulting in the formation of small clusters.

System:

We ran the tests on a machine with a Tesla K40c GPU and a Intel(R) Xeon(R) CPU E5-2640 v3 @ 2.60GHz. The tests for both the sequential and parallel implementations were made on the same machine.

Compilation:

We compiled the program using nvcc V7.5.17 with compiler options --std=c++11 -gencode arch=compute_30,code=sm_30. For the sequential code, we used g++ 4.9.2 with options --std=c++11 -O3.

Metric:

We ran simulations for 100 iterations, long enough to ensure that particle collisions happen frequently, except for the first iterations where the bodies are separated from each other. There we compute the average execution time and iteration averages for overlap correction cycles and edge flips.

The average execution time per time step, presented in Figure 9, shows two interesting features. First, for the case with long range forces the execution time is $O(N^{2})$, while for short range forces it is $O(N)$. Since both include the overlap correction algorithm, this implies that the execution time for the later is $O(N)$. Second, except for small systems, in the case of long range forces the execution time does not depend on the configuration, while for the short range forces, there is a clear dependence, with increasing complexity for $c_{4},c_{1},c_{2},c_{0},c_{3}$ (the same order of complexity is observed for long range forces at small $N$). This result is consistent with a cost $O(N)$ for the overlap correction, with a prefactor that may depend on the density and extension of the clusters.

To study the dependence of different configurations on the complexity of the overlap correction, in Figure 10 we plot the average per time step of iterations needed to correct all overlaps. The increasing complexity for $c_{4},c_{1},c_{2},c_{0},c_{3}$ is consistent with the previous results on Figure 9, because in $c_{3}$ most of the particles participate in corrections, while in $c_{2}$ almost half of the particles are excluded due to repulsion between same-type particles. However, this does not explain why $c_{4}$ has the least complexity factor, even though half of the particles overlap. This happens because all the overlaps on $c_{4}$ are corrected on the first iteration, which is probably due to the small size of the clusters. The number of iterations follow the same order in complexity as the execution time. Except for $c_{4}$ where clusters are disconnected, the number of iterations grow with $N$. This effect is due to the percolation of the large clusters, which cover the entire box and, therefore, the corrections become non-local and system size dependent. This growth is nevertheless weak, following an approximate logarithmic law. It is also noteworthy that the curves for long range forces are constant on $c_{4},c_{1}$ and $c_{0}$, slightly grow on $c_{2}$ and is relatively greater on $c_{3}$.

We also measured the performance of the Delaunay triangulation update algorithm, reporting the average of edge flip iterations made for both inverted triangle corrections and Lawson’s algorithm. The curves obtained on Figure 11 are less regular than the previous results, but keep the same general tendency. Unlike the curves for overlap correction, on where $c_{4}$ shows much smaller values than the other configurations, here the curve is comparable to $c_{0}$ and $c_{1}$. This happens because the underlying triangulation for $c_{4}$ has a great number of slivers, formed by the small particle density that forms relatively long edges. Then, according to the inverted triangle condition on section 4, it is more likely for $c_{4}$ to produce inverted triangle than the other configuration, whose triangles are more equilateral. The average for edge-flip iterations for long range forces has linear growth for all configurations, noting that $c_{2}$ and $c_{3}$ are the hardest cases to solve, as is the case on Figure 10. Though the number of iterations grows with $N$, it still remains negligible regarding the total time of a time step, so it is not a priority target for optimization.

Finally, we compared the n-body algorithms for long range forces, used on the different implementations without considering overlap corrections. shuffle is the presented optimization using warp-shuffle, while sharedMem is the GPU device memory algorithm described on section 4, observing a performance improvement of up to $2.4$ times from optimizing the parallel implementation for all tested values of $N$. It is also noteworthy that the optimized n-body algorithm allows simulation of $N=10^{6}$ particles at the same time that the sequential implementation solves the problem for $10^{5}$ bodies. For input sizes relevant to this study ($N\geq 10^{4}$), the time used by the sequential implementation is two orders of magnitude higher than the parallel solution, which allows the simulation of bigger particle systems for a longer physical time.

7 Validation

In order to verify that the developed overlap correction algorithm is efficient enough, we made two validation experiments. First, we test the locality of the correction, that is, how far it propagates though the system. For configurations $c0$ and $c3$ we print the result after $10^{5}$ time step iterations, painting with red the particles that took part in overlaps during the last simulated time step. Particles that did not take part in overlaps were painted green, so that every particle has a color. We repeat the process for decreasing values of $\Delta t$, expecting that the number of overlaps will decrease as the time step produces smaller movements. The results on Figure 13 allows us to verify the complexity factor associated to each configuration that shows up on the previous performance curves. Configuration $c0$ involves much less particles on overlap corrections than $c3$ upon lowering the time step. This decrease in execution time by reducing $\Delta t$ does not compensate, however, for the larger number of steps that are needed to achieve an specified physical time.

The second validation consists in testing the overlap correction on another colloidal model. We consider the Active Brownian Particle (ABP) model [4], where particles move in 2D with velocities of fixed magnitude $V_{0}$, with a direction that is specified by the director angle $\theta_{i}$. The integration rule for the positions after an interval $\Delta t$ is:

[TABLE]

In the same time interval, the angles $\theta_{i}$ are subjected to diffusive rotational Brownian motion, of amplitude $D$, and therefore evolve as:

[TABLE]

where $n_{i}$ is a random Gaussian variable of zero mean and unit variance.

We simulate the system with the same parameters used in Ref. [35], for two different packing fractions, obtaining the same phenomenlogy. At large packing fractions, the system evolves to the formation of a dense percolating dynamic cluster (see Fig. 14). Reducing the packing fraction, small clusters form, which merge in a slow coarsening process in the course of time as shown in Fig. 15.

8 Discussion

We presented algorithms for simulating colloidal particles subject to Brownian motion, interacting with short or long range force interactions, and presenting excluded volume. The overlap correction algorithm using Delaunay triangulations is a novel method. The algorithms implemented in CUDA for simulation are fully parallel, transferring data back to host only for measurements or outputs. The overlap correction algorithm can be used independently from the forces calculation, allowing to simulate different colloidal models including changed particles or self-propelled active systems. The Delaunay triangulation and the parallel edge-flip algorithm proved to be useful for solving overlaps efficiently. This opens the possibilty for using the Delaunay triangulation for solving related problems in the simulation, such as short range force calculation or approximated n-body simulations. The parallel n-body implementation was also successfully adapted and optimized to the particular conditions of colloidal particles, which opens up simulations of up to two orders of magnitude the number of particles used on the previous sequential implementation.

Acknowledgements

The authors would like to thank the NVIDIA GPU Research Center of the Department of Computer Science of the Universidad de Chile for supplying the equipment used for the tests presented here. This work was partially supported by the FONDECYT projects No. 1140778 and No. 3160182, and by project No. ENL009/15, VID, Universidad de Chile.