The Anderson metal-insulator transition (MIT) is central to our understanding of the quantum mechanical nature of disordered materials. Despite extensive efforts by theory and experiment, there is still no agreement on the value of the critical exponent ν describing the universality of the transition -the so-called "exponent puzzle". In this work, going beyond the standard Anderson model, we employ ab initio methods to study the MIT in a realistic model of a doped semiconductor. We use linear-scaling DFT to simulate prototypes of sulfur-doped silicon (Si:S). From these we build larger tight-binding models close to the critical concentration of the MIT. When the dopant concentration is increased, an impurity band forms and eventually delocalizes. We characterize the MIT via multifractal finite-size scaling, obtaining the phase diagram and estimates of ν. Our results suggest an explanation of the long-standing exponent puzzle, which we link to the hybridization of conduction and impurity bands. PACS numbers: 71.30.+h, The Anderson metal-insulator transition (MIT) is the paradigmatic quantum phase transition, resulting from spatial localization of the electronic wave function due to increasing disorder [1]. As for any such transition, universal critical exponents capture its underlying fundamental symmetries. This universality allows to disregard microscopic detail and the Anderson MIT is expected to share a single set of exponents. The last decade has witnessed many ground-breaking experiments designed to observe Anderson localization directly: with light [2-11], photonic crystals [9, 12], ultrasound [13, 14], matter waves [15], Bose-Einstein condensates [16] and ultracold matter [17,18]. The mobility edge [19], separating extended from localized states, was only measured for the first time in 2015 [20]. The hallmark of these experiments is the tunability of the experimental parameters and the ability to study systems where many-body interactions are absent or can be neglected. Under such controlled conditions, the observed exponential wave function decay, the existence of mobility edges and the critical properties of the transition [21,22] are in excellent agreement with the non-interacting Anderson model [1]. Furthermore, scaling at the transition [23] leads to high-precision estimates of the universal critical exponent ν from transport simulations (ν = 1.57(1.55, 1.59) [7]) and wave function statistics (ν = 1.590(1.579, 1.602) [4]). arXiv:1710.01742v4 [cond-mat.dis-nn]