We document a new method for including collisional return-to-isotropy in particle velocity distributions calculated by the MP-PIC method for numerical simulation of particle/fluid flows. Mathematically, we include the new method by adding to the transport equation for the particle distribution function (PDF), a Bhatnager, Gross, and Krook (BGK) collision term that causes velocity distributions to relax to isotropic Gaussian distributions. The numerical implementation is by a splitting technique in which we randomly sample from velocity distributions obtained as solutions to the PDF transport equation with just the return-to-isotropy source. Thus, collisions cause numerical particles to scatter in MP-PIC calculations, just as physical particles do in particle/fluid flows. The method is implemented in the Barracuda code, and two computational examples verify proper implementation of the method and show that more realistic results are obtained in calculations of impinging particle jets. In an appendix, we derive the particle continuum-flow (PCF) equations implied by the MP-PIC method in the collision dominated limit, and we obtain the collisional relaxation time for return-to-isotropy by matching the shear viscosity of the MP-PIC PCF equations, and the kinetic part of the shear viscosity used in PCF equations in the literature. (C) 2012 Elsevier Ltd. All rights reserved.