A new three-dimensional boundary-integral algorithm for deformable drops moving in a viscous medium at low Reynolds numbers is developed, which overcomes some familiar difficulties with boundary-integral calculations. The algorithm is used to simulate different modes of interaction between drops or bubbles, primarily for buoyancy-driven motion. The present iterative method for mean curvature calculation is found to be more robust and accurate than contour integration schemes. A novel iterative strategy based on combining biconjugate gradient and simple iterations overcomes the poor convergence of "successive substitutions" for drops in very close approach with extreme viscosity ratio. A substantially new variational method of global mesh stabilization solves the problem of mesh degradation with advantageous, soft stability constraints. A curvatureless boundary-integral formulation is also derived and shown to provide, in principle, a more accurate description of the drop breakup than the conventional formulation. The efficiency of these techniques is demonstrated by numerical examples for two drops in gravity-induced motion with high surface resolutions. The present code successfully simulates mutual approach of slightly deformable drops to extremely small separations, as well as their rotation when in "apparent contact," thus bridging the gap between finite deformation calculations and a recent asymptotic theory for small capillary numbers. Also provided is a 3D simulation of the experimental phenomenon of enhanced bubble coalescence, discovered by Manga and Stone [J. Fluid Mech. 256, 647 (1993); 300, 231 (1995)]. For drops of viscosity comparable to that of the surrounding fluid, it is shown in contrast that breakup is a typical result of hydrodynamic interaction in gravity-induced motion for large and even moderate capillary numbers. The code is readily applicable to any type of an ambient flow and may be adapted to more than two drops.