The importance of ceria () in many applications originates from the ease of oxygen vacancy formation and healing. The ordering of vacancies and the whereabouts of the excess charge in bulk are of no less significance than at ceria surfaces, but they have not received the same attention. In this work, the formation of neutral oxygen vacancies in bulk is investigated using density-functional theory (DFT) in the ( is an effective onsite Coulomb interaction parameter) approach for a broad range of vacancy concentrations . We find that the excess charge prefers to be localized in cation sites such that the mean coordination number is maximized, and if nearest-neighbor cation sites are reduced, they rather be nonuniformly distributed. Furthermore, we show that a vacancy repels other vacancies from its nearest-neighbor shell and that the [110] and [111] directions are possible directions for clustering of second- and third-neighbor vacancies, respectively. Vacancies prefer not to share cations. The results are discussed in a simple physical picture which enables the separation of the different contributions to the averaged vacancy formation energy. We also consider cells with fluorite structure and same stoichiometries as in existing bulk phases, i.e., (), (), and (), as well as the corresponding real structures. We find that the vacancy ordering and the location of the excess electrons are consistent with the results for single-phase reduced , but the , , and structures are substantially more stable. The stability of these phases as a function of pressure and temperature is discussed. Vacancy-induced lattice relaxation effects are crucial for the interpretation of the results.