Incorporating a link-cell search into the velocity-Verlet algorithm gives, in outline,
Read in initial data.
Read in configurational data (solid inclusions).
Set up the map to find neighboring cells (subroutine MAPS).
Perform the QDPD cycle for each time step.
Given the forces ƒ1i acting on particles at time t, the fundamental QDPD cycle, repeated for as many timesteps as are in a simulation, is
Compute the new particle positions from
Compute the midpoint velocity (velocity at the midpoint of the time step) from
Create the linked list (subroutines TOPMAP and LINKS)
Calculate new forces ƒ2i.
Compute new velocities from the new forces
And then the cycle begins again.
In the basic QDPD (MD) cycle above, the r2s may belong to particles which have moved out of the simulation box, such as particles in the periodic image of cell 21 represented by the dashed box 21' in Fig. 2. This can be handled by introducing another set of coordinates, r3, given by
where Lx, Ly, and Lz are the simulation dimensions and ANINT is the function which returns the nearest whole integer to its argument. The r3 coordinates are used in creating linked lists of particles in cells prior to the force calculations. Consequently, particles which have moved into 21' will end up assigned to 21 which is where the formula for nearest neighbors expects to find them. Hence the r3 coordinates make sure that particles are assigned to one of the cells in the QDPD simulation box (i.e., particles are kept within the QDPD simulation box running from (−Lx/2, −Ly /2, −Lz /2) to (Lx /2, Ly /2, Lz /2) in r3 space). They also are consistent with the proper mapping of nearest neighbor cells given by the ICELL(Ix, Iy) formula. One other point has to do with calculating forces on particles > rc away (since the particles may have moved out of the simulation box). In calculating forces, r2 coordinates are used, and these r2 coordinates may be > r2 away, a violation of the minimum image condition. To correct for this, the difference between particles in the forces calculation is
Consider Fig. 2 again. With the ANINT corrections above, particles in the cells 4, 5, 1, and 21 are within r2 away from particles in cell 25. This is the way our sequential version of the program was written. Our parallel version of the program does this differently. In treating edge cells, a "ghost" layer of cells is added to the QDPD simulation box. The dashed cells in Fig. 2, the nearest neighbor periodic cells, are part of the "ghost" layer of cells. The formation of these "ghost" cells will be discussed in the next section.
In the sequential version of the program, MAPS considers all cells in all layers except for the top layer (the topmost layer in Y in a 3D simulation), and computes only half of these cells, taking into account Newton's third law. Because QDPD in its present form is being used to study the steady-shear viscosity of a suspension of solid inclusions in a Newtonian fluid, there is a shear boundary condition at the topmost layer of the QDPD simulation box, implemented with the Lees-Edwards boundary conditions [9] (pp. 246-247). These boundary conditions simulate a uniform shear in the XY plane (i.e, a constant velocity gradient is set up in the Y direction and the actual shear occurs in the X direction). Figure 3 shows a time series of the motion of a single ellipsoidal inclusion subject to shear. Proceeding from left to right, the different colors (or greyscale levels) [12] correspond to the time sequence. The single ellipsoid rotation is a well known phenomenon seen in experiments called Jeffery's orbits. The shearing boundary conditions were obtained by applying a constant strain rate to the right at the top of the figure and to the left at the bottom of the figure. Figure 4 [9] (p. 246) demonstrates this situation in the context of the computer simulation. The central box in the figure is the QDPD simulation box (the entire box in Fig. 5, not just the one on the central processor 4). Boxes in the layer above are moving at a certain speed in the positive direction, and boxes in the layer below are moving at the same speed in the negative direction. To implement this shear boundary condition at the topmost layer (because of Newton's third law, we only have to treat the top), the top layer is tackled separately in subroutine TOPMAP in our sequential version of the program. Its purpose is to create the list of neighboring cells for the topmost layer of cells, taking into account the movement of the cells with respect to each other due to shear. TOPMAP is called every timestep in the simulation just before the force calculation, but only on the topmost processors. The periodic minimum image convention must also be modified to account for this shear. The r3s are modified to be where the upper layer (BCD in Fig. 4) is displaced relative to the central box by an amount strain. Similar corrections are made in forces.
Fig. 3. Tumbling of a single ellipsoidal inclusion under shear. Details of the algorithm for a rigid body are given in [5].
Fig. 4. Lees-Edwards boundary conditions for homogeneous shear (adopted from Allen and Tildesley, Computer Simulation of Liquids, Oxford, 1987, Fig. 8.2).
Fig. 5. A 9 processor 2-D domain. The small rectangles are cells associated with the link-cell algorithm. The dashed lines correspond to the ghost cells.
