A method for the study of positronium annihilation in polymers is presented. The polymer configurations are generated by atomistic molecular dynamics simulations; their geometry is fully known, so that the size and the shape of the holes can be analyzed. The positronium lifetime in the free volume sites is calculated at finite temperature with the path integral Monte Carlo method in the framework of an extended model, which comprises the positronium-polymer interaction and the electron density in the polymer. The spectra and free volume properties are compared to results from positronium annihilation spectroscopy experiments. As an example, the method is applied to polystyrene at ambient conditions.