Observations of surface emission from isolated neutron stars (NSs) provide unique challenges to theoretical modelling of radiative transfer in magnetized NS atmospheres. Recent work has demonstrated the critical role of vacuum polarization effects in determining NS spectra and polarization signals, in particular the conversion of photon modes (due to the 'vacuum resonance' between plasma and vacuum polarization) propagating in the density gradient of the NS atmosphere. Previous NS atmosphere models incorporated the mode conversion effect approximately, relying on transfer equations for the photon modes. Such treatments are inadequate near the vacuum resonance, particularly for magnetic field strengths around B ∼ B l ≃ 7 × 10 13 G, where the vacuum resonance occurs near the photosphere. In this paper, we provide an accurate treatment of the mode conversion effect in magnetized NS atmosphere models, employing both the modal radiative transfer equations coupled with an accurate mode conversion probability at the vacuum resonance, and the full evolution equations for the photon Stokes parameters. In doing so, we are able to quantitatively calculate the effects of vacuum polarization on atmosphere structure, emission spectra and beam patterns, and polarizations for the entire range of magnetic field strengths, B = 10 12 -10 15 G. In agreement with previous works, we find that for NSs with magnetic field strength B ≥ 2B l , vacuum polarization reduces the widths of spectral features, and softens the hard spectral tail typical of magnetized atmosphere models. For B ≤ B l /2, vacuum polarization does not change the emission spectra, but can significantly affect the polarization signals. Our new, accurate treatment of vacuum polarization is particularly important for quantitative modelling of NS atmospheres with 'intermediate' magnetic fields, B ≃ (0.5-2) B l . We provide fitting formulae for the temperature profiles for a suite of NS atmosphere models with different field strengths, effective temperatures and chemical compositions (ionized H or He). These analytical profiles are useful for direct modelling of various observed properties of NS surface emission. As an example, we calculate the observed intensity and polarization light curves from a rotating NS hotspot, taking into account the evolution of photon polarization in the magnetosphere. We show that vacuum polarization induces a unique energy-dependent linear polarization signature, and that circular polarization can be generated in the magnetosphere of rapidly rotating NSs. We discuss the implications of our results to recent observations of thermally emitting isolated NSs and magnetars, as well as the prospects of future spectral and polarization observations.