A simple and efficient internal-coordinate importance sampling protocol for the Monte Carlo computation of (up to fourth-order) virial coefficients ̅B(n) of atomic systems is proposed. The key feature is a multivariate sampling distribution that mimics the product structure of the dominating pairwise-additive parts of the ̅B(n). This scheme is shown to be competitive over routine numerical methods and, as a proof of principle, applied to neon: The second, third, and fourth virial coefficients of neon as well as equation-of-state data are computed from ab initio two- and three-body potentials; four-body contributions are found to be insignificant. Kirkwood-Wigner quantum corrections to first order are found to be crucial to the observed agreement with recent ab initio and experimental reference data sets but are likely inadequate at very low temperatures.