What are the advantages and disadvantages of the particle decomposition and domain decomposition parallelization algorithms?

Question

I am running molecular dynamics (MD) simulations using several software packages, like Gromacs and DL_POLY.

Gromacs now supports both the particle decomposition and domain decomposition algorithms. By default, Gromacs simulations use domain decomposition, although for many years, until recently, particle decomposition was the only method implemented in Gromacs. In one of the Gromacs papers (DOI 10.1002/jcc.20291), the authors give a reason for their initial choice of particle decomposition:

"An early design decision was the choice to work with particle decomposition rather than domain decomposition to distribute work over the processors. In the latter case, spatial domains are assigned to processors, which enables finding spatial neighbors quickly by local communication only, but complications due to particles that move over spatial boundaries are considerable. Domain decomposition is a better choice only when linear system size considerably exceeds the range of interaction, which is seldom the case in molecular dynamics. With particle decomposition each processor computes the forces and coordinate/velocity updates for an assigned fraction of the particles, using a precomputed neighborlist evenly distributed over processors. The force $F_{ij}$ arising from the pair interaction between particles $i$ and $j$, which is needed for the velocity update of both particles $i$ and $j$, is computed only once and communicated to other processors. Every processor keeps in its local memory the complete coordinate set of the system rather than restricting storage to the coordinates it needs. This is simpler and saves communication overhead, while the memory claim is usually not a limiting factor at all, even for millions of particles. The neighborlist, on the other hand, which can contain up to 1000 times the number of particles, is distributed over the processors. Communication is essentially restricted to sending coordinates and forces once per time step around the processor ring. These choices have proven to be robust over time and easily applicable to modern processor clusters."

What do they mean by "linear system size" in the sentence "Domain decomposition is a better choice only when linear system size considerably exceeds the range of interaction, which is seldom the case in molecular dynamics"? From the paragraph above, I get the idea that particle decomposition has the advantage that one does not have to deal with particles moving across domain boundaries; rather, you just have to have enough memory for each processor to store the total system configuration. So particle decomposition is looking very favorable, whereas domain decomposition is looking very unfavorable.

I am sure that this is a very complicated question (and probably the subject of many books), but just basically, if particle decomposition seems so favorable, why would anyone need to use domain decomposition? Is domain decomposition just favorable if the system's size is very large (making it difficult or impossible to store the total configuration in each processor)? Based on the quoted paragraph above, I am not sure why domain decomposition is now, just recently, the default parallelization algorithm in Gromacs.

It seems that DL_POLY now (version 4) also uses domain decomposition. From the version 4 manual:

"The division of the conguration data in this way is based on the location of the atoms in the simulation cell, such a geometric allocation of system data is the hallmark of DD algorithms. Note that in order for this strategy to work efficiently, the simulated system must possess a reasonably uniform density, so that each processor is allocated almost an equal portion of atom data (as much as possible). Through this approach the forces computation and integration of the equations of motion are shared (reasonably) equally between processors and to a large extent can be computed independently on each processor. The method is conceptually simple though tricky to program and is particularly suited to large scale simulations, where efficiency is highest.

...

In the case of the DD strategy the SHAKE (RATTLE) algorithm is simpler than for the Replicated Data method of DL_POLY Classic), where global updates of the atom positions (merging and splicing) are required."

This makes it sound as if domain decomposition is good because it may be more efficient, even though perhaps more difficult to implement.

On the other hand, a previous version (DL_POLY Classic) used replicated data parallelization, which seems to be another name for particle decomposition. From that version's manual:

The Replicated Data (RD) strategy is one of several ways to achieve parallelisation in MD. Its name derives from the replication of the configuration data on each node of a parallel computer (i.e. the arrays defining the atomic coordinates $\textbf{r}_i$, velocities $\textbf{v}_i$, and forces $\textbf{f}_i$, for all $N$ atoms in the simulated system, are reproduced on every processing node). In this strategy most of the forces computation and integration of the equations of motion can be shared easily and equally between nodes and to a large extent be processed independently on each node. The method is relatively simple to program and is reasonably efficient. Moreover, it can be “collapsed” to run on a single processor very easily. However the strategy can be expensive in memory and have high communication overheads, but overall it has proven to be successful over a wide range of applications.

This paragraph seems generally consistant with the first paragraph in this question, except that it says that replicated data/particle decomposition has "high communication overheads." The paragraph from the Gromacs paper seems to say the opposite -- that particle decomposition is preferable because it has lower communication overhead than domain decomposition.

Do you have any thoughts?

Answer 1

Particle and domain decomposition are directly connected to the two main methods of speeding up force calculations for systems with limited-range interactions - Verlet neighbour lists and cell linked lists. If you'd like to get into details, there is a pretty nice book from Allen and Tildesley, called Computer Simulation of Liquids, considered by many to be the "bible" of Molecular Dynamics and Monte Carlo studies. Then there is the Numerical Simulation in Molecular Dynamics from Griebel, Knapek and Zumbusch, that goes deep into the various techniques for parallel implementation of MD.

Basically Verlet lists build a list of all the neighbours of a given atom/molecule (or particles in general) within a given radius. Later when pairs of atoms are being examined in order to compute the force, the list is consulted. Once that you have the list constructed it is obvious which particles are close to which other and they can be distributed among different processors for evaluation. The list is constructed only every now and then with some clever techniques to keep it up to date, since it takes $O(N^2)$ to build (all possible paris of particles are examined). Force evaluation given the already constructed list takes $O(N)$.

Cell linked lists divide the space in equally sized cells that are larger than the cut-off distance of the interaction potential and then each particle is put on a list associated with the cell it falls into. This process takes $O(N)$. Then neighbours of a given particle are only searched for in the same cell or in its neighbouring cells. Since each cell has a constant number of neighbours (e.g. 26 in the 3-D case), forces are evaluated in $O(N)$. But here the constant multiplier could be large enough so to make this algorithm to scale worse than the Verlet list method. But for large enough $N$ it scales better. Hence the linear size argument. The domain decomposition method is a straightforward extension of the cell linked lists method - cells are divided among different CPUs.

The problem with domain decomposition is that it has to communicate when particles move from one cell to another one that is taken care of by another CPU. This could become problematic at higher simulation temperatures where particles tend to move further away than their equilibrium position, or when there is a flow of particles. Also information from cells on the domain border has to be transfered on each iteration to the neighbouring domain(s). But all this is locally synchronous communication and could be done very efficiently.

The replicated data is the easiest approach but unfortunately it requires that at every step all position and velocity information is synced globally. This really doesn't scale well and for very large system the global ammount of memory is the size of the data structure times the number of CPUs used, while one of the goals of the parallel processing is distribution of data such that each CPU holds less than the global ammount of data.

In summary, there exists no "one size fits all" method, suitable for all systems being simulated. Most of the time the best parallelisation strategy can be deduced from the system geometry and the appropriate for that case MD code could be picked - they all implement more of less the same underlying force fields and integrators after all.

Answer 2

By "Domain decomposition is a better choice only when linear system size considerably exceeds the range of interaction, which is seldom the case in molecular dynamics" the authors of that (very old) GROMACS paper mean that if the spatial size of the neighbour list is of the order of 1 nm, and the simulation cell is only several nanometers, then the overhead from doing domain decomposition is too high. You may as well accept an all-to-all information distribution in particle decomposition, and not need to spend time on all the book-keeping for domain decomposition.

The problem with particle decomposition as GROMACS implemented it was that over time the particles assigned to each processor diffuse through space. Since responsibility for computing each interaction was fixed by their initial location, the diffusion gradually increased the volume of the total space each processor needed to know in order to build its neighbour list, even if the total computation described by the neighbour list was constant. In practice, you would periodically re-start the simulation to reset the data and communication locality.

Your supposition that "particle decomposition has the advantage that one does not have to deal with particles moving across domain boundaries" does not hold if diffusion is significant over the time scale of the simulation.

Domain decomposition deals with this "up front" by migrating responsibility for the interaction along with the diffusion, thereby improving data locality on each processor, and minimizing communication volume.

Disclaimer: I help develop GROMACS, and will probably rip out the particle decomposition implementation next week ;-)

Answer 3

I would like to add to the answer of Hristo Iliev. While his post mostly talks about the computational complexity, when it comes to parallelization, the communication complexity is at least as important - and that it the main reason for domain decomposition.

Modern parallel machines usually have some kind of a torus topology. This means that each CPU has a number of "neighboring" CPUs that it can communicate to very quickly. Communicating to a CPU that is not a neighbor is more costly. Therefore, it is always favorable to have an algorithm that only needs to communicate with neighboring CPUs.

When using a particle decomposition, the interaction partners of a particle are randomly distributed on all other CPUs. To be able to compute the interactions, it needs to know the coordinates of all partners, so it needs to communicate with all other CPUs. Effectively, this means that when you have $P$ CPUs, $\mathcal{O}(P^2)$ communication steps are required. Furthermore, a lot of the communication is with non-neighboring CPUs.

On the other hand, in a domain decomposition scheme, all interaction partners live on a neighboring CPU. This means that it needs to communicate only with its nearest neighbors to get the updated information, i.e. $P$ CPUs require $\mathcal{O}(P)$ communication steps. Furthermore, all of the communication is with neighboring CPUs.

Clearly, when the system is non-uniformly distributed, this scheme doesn't work as optimal. When that happens, this either means that some CPUs have significantly more work than others, or that one has to dynamically adapt the domains. This might cause a mismatch between the domain topology and the network topology (large domains have more neighbors). Still, the communication complexity of $\mathcal{O}(P)$ still holds.

Note, however, that non-uniform systems are not as common as it may sound, they only occur when either simulating something in vacuum, or when using an implicit solvent. The densities of crystals and liquids are close enough to run domain decomposition.