We investigate the dynamics of inertial particles in homogeneous isotropic turbulence, under one-way momentum coupling, using a new computational approach that incorporates the effect of long-range many-body aerodynamic interactions along with the short-range lubrication forces. The implementation couples hybrid direct numerical simulations (HDNS) with the analytical solutions of two rigid spheres moving in an unbounded fluid. Concerning the velocity field seen by the particles, the algorithm switches from the flow solution in terms of HDNS to analytical formulae when the separation distance between particles becomes comparable to their average radius. Standard HDNS is unable to correctly represent the short-range interactions since this method is based on the superposition of the Stokes solutions for single spheres. Our results show that for the turbulent kinetic energy dissipation rates typical of atmospheric clouds, the radial relative velocities (RRVs) of the droplets increase, and the radial distribution function (RDF) decreases in the near-contact region if the lubrication forces are taken into account. These changes are more pronounced when the effect of gravity is considered. Away from the contact region, there is not much change in RRVs and RDFs. For turbulent clouds with lower dissipation rates lubrication forces significantly enhance the average RRV in the limit of low Stokes number. This enhancement, however, is statistically insignificant because the number of particle pairs at close proximity is very small. The effect of mass loading on the collision statistics is also investigated, demonstrating an increase in RRV and a reduction in RDF with the droplet concentration.