A Lagrangian method of tracking particles in a fluid medium is conceptually simple and numerically robust. It is more practical and efficient than two-fluid methods for dilute suspensions of particulate phase, where particle-fluid and particle-particle interactions (four-way coupling) are weak [{Elghobashi, }{1994}]. However, when it comes to dense suspensions it stumbles on technical problems of computing particle interactions. For each particle, identification of neighboring particles and expensive looping over them become necessary. This may considerably degrade the efficiency of the Lagrangian method, and two fluid methods are usually favored in this case. Nevertheless, in many practical situations, like spray and coal combustion, dynamics of aerosols and particulate pollutants, bubbly flows, etc. concentration of the particulate phase is very uneven with small regions of very dense suspension and large voids with very little or no particles. Lagrangian particle tracking algorithms are well suited to handle this unevenness, offering savings in computer memory and execution time. Therefore, it is important to develop an efficient four-way coupling technique that would account for particle-fluid and particle-particle interactions in a Lagrangian particle tracking method.