There has been recent interest in mesoscopic models of hydrodynamics called dissipative particle dynamics (DPD) that blend cellular automata ideas with molecular dynamics methods [1]. The original DPD alogrithm utilized symmetry properties such as conservation of mass, momentum and Galilean invariance to construct a set of equations for updating the position of particles which can be thought of as representing clusters of molecules or "lumps" of fluid. Later modifications of the DPD algorithm resulted in a more rigourous formulation which was consistent with the fluctuation dissipation theorem [2]. Improvements to the temperature behavior of the DPD algorithm were made by incorporating a velocity Verlet algorithm which allowed a larger time step while still producing a satisfactory temperature control [3]. An algorithm for modeling the motion of arbitrary shaped objects subject to hydrodynamic interactions by DPD was suggested by Koelman and Hoogerbrugge [4]. The rigid body is approximated by "freezing" a set of randomly placed particles where the solid inclusion is located and updating their position according to the Euler equations. The original DPD alogrithm used an Euler algorithm for updating the postions of free particles and a leap frog algorithm for updating the postion of the rigid body. A motivation of this work was to develop an efficient algorithm to update both the free particles and the rigid body position in a manner consistent with the velocity Verlet algorithm.
A commonly used velocity Verlet based algorithm for updating the position of rigid bodies is the so called rattle algorithm [5]. The rattle routine solves the constraint equations of the particles comprising the rigid body by a relaxation method. Further, the stress tensor (from an atomic view) can be directly obtained from the constraint forces calculated in the algorithm and is completely symmetric. While the rattle routine is of order N 2 (where N is the number of particles in the inclusion) and, hence, can be prohibitively slow for the case of modeling the motion of solid inclusions composed of large numbers of particles, it serves as an accurate benchmark to test other algorithms.
In this brief report we show how the use of quaternions to represent the orientation of such objects [6,7] can greatly increase the computational efficiency of DPD simulations. Quaternions provide a convenient way to represent the orientation of rigid objects since the transformations between body-fixed and laboratory coordinate reference frames contain no singularities when expressed as quaternions. First, we review the development of the equations of motion for the quaternions. Next, we indicate how to efficiently apply the velocity Verlet algorithm [8,9] to the quaternion equations [10] and demonstrate its use in the simulation of water. We then discuss the modifications of the algorithm needed to include the velocity dependent dissipative forces in DPD simulations. A procedure for determining the rigid body's contribution to the stress tensor, consistent with the velocity Verlet algorithm, is given and compared to that derived from the rattle routine.