DPD motion of large rigid bodies

Forces in DPD depend both on the relative positions and velocities of particles. Hence a predictor velocity must be estimated to input into the force calculations. A reasonable approach to determine a predictor velocity could be based on an estimate of angular velocity at t = $t=\delta t/2$ t/2 derived from the quaternion equations of motion. However, to more closely match the trajectories obtained from the rattle routine that strictly follow the velocity-Verlet algorithm, the predictor velocity was simply based on the average velocity obtained in moving from postion x(0) to x( $x(\delta t)$ t ) in Eq. (9). Otherwise, the quaternion and rattle routines may not be consistent except in the limit of infinitesimally small time step.

The following modification was then made to velocity equation in Eq. (9):

$$\tilde{\dot x}({\delta t\over 2}) = \dot x(0) + {\lambda \, \delta t}\, a(0) \approx {x(\delta t)- x(0) \over \delta t },$$ (15)

and

$$\dot x(\delta t) = \dot x(0) + {\delta t\over 2}[ a(0) + a (\delta t,\tilde{\dot x} ({\delta t\over 2}))].$$ (16)

for the rattle routine and for the center of mass motion of the rigid body when using the quaternion-based algorithm. The final velocity of the constituent particles of the rigid body was derived from the angular velocity at the end of the time step. We used $\lambda={1\over2}$ = ½ for our simulations. For further discussion of the effect of varying $\lambda$ , see Ref.[3].

While the rigid body contribution to the stress tensor is readily calculated from the rattle routine, the constraint forces contributions are not immediately obtained from our quaternion algorithm. However, because the motion of the rigid body closely follows the trajectory obtained by the rattle routine, we can approximate the constraint forces contributions by considering the velocity Verlet algorithm. Let a_i (0) be the acceleration of particle i on the rigid body which results from the sum of non-constraint forces due to all particles (including those in the rigid body) and the constraint forces from particles in the rigid body. That is, $m\vec{a}_i =\vec{F}_i={\vec{F}^{nc}}_i+{\vec{F}^c}_i$ where the superscripts nc and c correspond to non-constraint forces and constraint forces, respectively. Since all the non-constraint forces are known, as well as the velocity and positions of the particle at t = 0 and t = $t=\delta t$ t, the sum of the constraint forces on particle i can be derived by assuming that particles follow the same time evolution as that derived from the rattle routine using Eqs. (9) and (16).

We now show that the sum of the constraint forces on each particle is all that is needed to correctly determine the constraint force contribution to the stress tensor. First, the contribution to the stress tensor from constraint forces is:

$$\sigma_{\alpha \beta}^c = {1\over 2} \sum_{i,j} {F_{ij}^c}_\alpha(\vec{r}_i-\vec{r}_j)_\beta$$ (17)

where i and j correspond to all particles in the rigid body. Note that ${\vec{F}^c}_{ij}=-{\vec{F}^c}_{ji}$ and that the sum of constraint forces on particle i due to particles j is ${\vec{F}}_i^c=\sum_j {\vec{F}}_{ij}^c$. We can then write

$$\begin{eqnarray*}\sigma_{\alpha \beta}^c &= {1 \over 2} \sum_{i,j} {{{F}}_{ij}... ..._j}_\beta \\ &=\sum_{i} {{{F}}_{i}^c}_\alpha{\vec{r}_i}_\beta \end{eqnarray*}$$

Since $\vec{F}^{nc}$is known and $\vec{a}_i$ can be estimated from Eqs. (9) and (16), we can determine $\vec{F}^c$to the accuracy of the algorithm and determine contributions to the stress tensor from the two parts of the velocity Verlet algorithm. Comparing stress tensor values between the rattle and our quaternion algorithm, we found agreement to six significant figures and, further, that the stress tensor was symmetric to the same order. Indeed for the time step used here, $\delta t = 0.01$ t = 0.01, out results are consistent with ${\cal O}(\delta t^3)$, the accuracy of the velocity-Verlet algorithm.