The lattice Monte Carlo method with parameters from recent first-principle calculations1,2 are used to investigate dopant diffusion in silicon. In the simulations, vacancy hopping on a silicon lattice is biased by changes in system energy, including interactions up to the sixth-nearest neighbor. We find that vacancy-mediated diffusivity increases dramatically above 1020 cm−3, in agreement with experimental observations3 and previous calculations.4 However, for very long simulation times, arsenic diffusivity is reduced due to formation of AsxV complexes, with clustering more pronounced at high doping levels. As suggested by Ramamoorthy and Pantelides,5 we find that As2V complexes are mobile, and although they diffuse much more slowly than AsV pairs, they appear likely to have a significant role in high concentration diffusion due to their much higher numbers. We also investigated dopant fluxes in a vacancy gradient. For dopants like As for which pair diffusion is limited by the dissociation to third-nearest neighbor distances, the dopant flux is less than that predicted by pair diffusion models, with greater difference at higher temperatures. In contrast, for phosphorus/vacancy pairs, whose diffusion is limited by dopant/vacancy exchange, the dopant flux is close to the predictions of pair diffusion.