We present a simulation of the cosmic evolution of the atomic and molecular phases of the cold hydrogen gas in about 3x10^7 galaxies, obtained by postprocessing the virtual galaxy catalog produced by De Lucia & Blaizot on the Millennium Simulation of cosmic structure. Our method uses a set of physical prescriptions to assign neutral atomic hydrogen (HI) and molecular hydrogen (H2) to galaxies, based on their total cold gas masses and a few additional galaxy properties. These prescriptions are specially designed for large cosmological simulations, where, given current computational limitations, individual galaxies can only be represented by simplistic model objects with a few global properties. Our recipes allow us to (1) split total cold gas masses between HI, H2, and helium, (2) assign realistic sizes to both the HI and H2 disks, and (3) evaluate the corresponding velocity profiles and shapes of the characteristic radio emission lines. The results presented in this paper include the local HI and H2 mass functions, the CO luminosity function, the cold gas mass-diameter relation, and the Tully-Fisher relation (TFR), which all match recent observational data from the local universe. We also present high-redshift predictions of cold gas diameters and the TFR, both of which appear to evolve markedly with redshift.