A three-dimensional boundary-integral algorithm has been developed to handle the tangential Marangoni stresses in thermocapillary motion of drops. Depending on whether the integration and observation points are on the same or different drops, singularity or near-singularity subtraction is used in the inhomogeneous term of the boundary-integral formulation. Integration is then performed analytically over flat triangles in the subtracted region. Relative trajectories for two deformable drops are calculated for different values of the drop size ratio, drop-to-medium thermal conductivity ratio, and viscosity ratio and compared to those for spherical and slightly deformable drops. Results indicate that deformation increases the minimum separation and inhibits coalescence but is not important enough for appropriate physical parameters to induce the capture or breakup behaviors observed in buoyancy. Interaction times calculated by artificially continuing spherical drop trajectories yield results accurate to within about 10%.