We examined the underlying mechanisms for thermal conductivity suppression in crystalline silicon by substitutional doping with different elements (X = boron, aluminum, phosphorus, and arsenic). In particular, the relative effects of doping-induced mass disorder, bond disorder, and lattice strain were assessed using nonequilibrium molecular dynamics simulations. Stillinger-Weber potential parameters for SiX interatomic interactions were optimized by fitting to relevant atomic forces from first-principles calculations. We first calculated the thermal conductivity variation of B-doped Si as a function of dopant concentration; the result shows excellent agreement with existing experimental data, indicating the reliability of our force-field-based simulations. At the dopant concentration of about 5 × 10 20 cm −3 , the Si thermal conductivity value is predicted to be reduced from 137 W/mK at 300 K in undoped Si to 18/39/57/78 W/mK in As/B/P/Al-doped Si. Our study demonstrates that the mass disorder effect is primarily responsible for the thermal conductivity suppression in the As-and B-doped cases, whereas the bond disorder contribution is found to be more important than the mass disorder contribution in the Al-and P-doped cases; for all these systems, the lattice strain effect turns out to play a minor role in the reduction of lattice thermal conductivity.