We report an improved method for the calculation of tunneling splittings between degenerate configurations in molecules and clusters using path-integral molecular dynamics (PIMD). Starting from an expression involving a ratio of thermodynamic density matrices at the bottom of the symmetric wells, we use thermodynamic integration with molecular dynamics simulations and a Langevin thermostat to compute the splittings stochastically. The thermodynamic integration is performed by sampling along the semiclassical instanton path, which provides an efficient reaction coordinate as well as being physically well-motivated. This approach allows us to carry out PIMD calculations of the multi-well tunnelling splitting pattern in water dimer, and to refine previous PIMD calculations for one-dimensional models and malonaldehyde. The large (acceptor) splitting in water dimer agrees to within 20% of benchmark variational results, and the smaller splittings are within 10%.