We construct hydrogen atmosphere models for magnetized neutron stars in radiative equilibrium with surface fields B = 10 12 −5×10 14 G and effective temperatures T eff ∼ a few×10 5 −10 6 K by solving the full radiative transfer equations for both polarization modes in the magnetized hydrogen plasma. The atmospheres directly determine the characteristics of thermal emission from isolated neutron stars. We study the effects of vacuum polarization and bound atoms on the atmosphere structure and spectra. For the lower magnetic field models (B ∼ 10 12 G), the spectral features due to neutral atoms lie at extreme UV and very soft X-ray energies and therefore are not likely to be observed. However, the continuum flux is also different from the fully ionized case, especially at lower energies. For the higher magnetic field models, we find that vacuum polarization softens the high energy tail of the thermal spectrum. We show that this depression of continuum flux strongly suppresses not only the proton cyclotron line but also spectral features due to bound species; therefore spectral lines or features in thermal radiation are more difficult to observe when the neutron star magnetic field is > ∼ 10 14 G.