Positron annihilation in solid-state matter can be utilized to detect and identify open-volume defects. The momentum distribution of the annihilation radiation is an important observable in positron-based measurements and can reveal information on the chemical surroundings of the defect sites. In this paper, we present a variational quantum Monte Carlo method for simulation of the momentum densities of annihilating electron-positron pairs in semiconductors and insulators. We study finite-size effects, effects of lattice vibrations, and different levels of trial wave functions. Small simulation cells and simple wave function forms are found to be sufficient for accurate calculations in the simulation of pristine lattices, enabling cheap accumulation of results. We compare calculated predictions of the Doppler broadening of the 511-keV 2γ annihilation line in aluminium nitride and silicon against experimental data measured from reference samples. Our results achieve better agreement with experiments in these materials than conventional state-of-the-art methods and prove that direct modeling of the electron-positron correlations is important for a supporting theory of positron annihilation spectroscopies.