We present a new efficient calculation to propagate cosmic ray muons from the surface of the Earth to deep underground laboratories, allowing us to look at the physics and performance of various models of high-energy cosmic rays. The evolution of cosmic rays in the Earth's atmosphere is computed with MCEq (Matrix Cascade Equation), taking into account different combinations of primary and hadronic interaction models in order to calculate the muon flux at the surface. The latter serves as an input for the Monte Carlo code PROPOSAL (Propagator with Optimal Precision and Optimised Speed for All Leptons) to propagate the muons through the rock. A forward prediction for underground muon spectra at different slant depths, including the muon survival probabilities and underground energy spectra, is calculated with very high precision. The reliability of this state-of-the-art calculation was achieved by comparing the results obtained for the vertical muon intensity and total muon flux with the measured data at various underground sites with both flat overburdens and mountains. The implications of the results as well as the seasonal variation of the muon flux will also be discussed.