Current and upcoming radio telescopes will map the spatial distribution of cosmic neutral hydrogen (H I) through its 21 cm emission. In order to extract the maximum information from these surveys, accurate theoretical predictions are needed. We study the abundance and clustering properties of H I at redshifts z5 using TNG100, a large state-of-the-art magnetohydrodynamic simulation of a 75h −1 Mpc box size, which is part of the IllustrisTNG Project. We show that most of the H I lies within dark matter halos, and we provide fits for the halo H I mass function, i.e.,the mean H I mass hosted by a halo of mass M at redshift z. We find that only halos with circular velocities larger than ;30km s −1 contain H I. While the density profiles of H I exhibit a large halo-to-halo scatter, the mean profiles are universal across mass and redshift. The H I in low-mass halos is mostly located in the central galaxy, while in massive halos the H I is concentrated in the satellites. Our simulation reproduces the bias value of damped Lyα systems from observations. We show that the H I and matter density probability distribution functions differ significantly. Our results point out that for small halos, the H I bulk velocity goes in the same direction and has the same magnitude as the halo peculiar velocity, while in large halos, differences show up. We find that halo H I velocity dispersion follows a power law with halo mass. We find a complicated H I bias, with H I already becoming nonlinear at k=0.3 h Mpc −1 at z3. The clustering of H I can, however, be accurately reproduced by perturbative methods. We find a new secondary bias by showing that the clustering of halos depends not only on mass but also on H I content. We compute the amplitude of the H I shot noise and find that it is small at all redshifts, verifying the robustness of BAO measurements with 21 cm intensity mapping. We study the clustering of H I in redshift space and show that linear theory can explain the ratio between the monopoles in redshift and real space down to 0.3, 0.5, and 1 h Mpc −1 at redshifts 3, 4, and 5, respectively. We find that the amplitude of the Fingers-of-God effect is larger for H I than for matter, since H I is found only in halos above a certain mass. We point out that 21 cm maps can be created from N-body simulations rather than full hydrodynamic simulations. Modeling the one-halo term is crucial for achieving percent accuracy with respect to a full hydrodynamic treatment. Although our results are not converged against resolution, they are, however, very useful as we work at the resolution where the model parameters have been calibrated to reproduce galaxy properties.