Heteroaromatics and benzene derivatives constitute important structural classes in drug discovery, agrochemistry, and material science. Their halogenated derivatives are often applied as substrates in carbon-carbon and carbon-heteroatom cross-coupling reactions such as Suzuki-Miyaura, Heck, and Buchwald-Hartwig couplings. The prerequisite (hetero)aryl halides are typically prepared by electrophilic aromatic substitution (EAS). Halogenated arenes are also important substrates for directed metallation reactions and for the generation of organolithium and organomagnesium species by metal-halogen exchange, metal insertion, or direct metallation. Unlike the benzene series for which the relative reactivity of the substrates is textbook knowledge, it is not a priori obvious at which position(s) halogenation will occur for heteroaromatic systems, especially in compounds that contain multiple (hetero)aromatic rings or in compounds that contain both heteroarene and benzene rings. Consequently, organic chemists tend to install the halogens early in the synthesis because they are not comfortable with late-stage functionalization. Thus, there is an unmet need for synthetic chemists to be able to predict the regioselectivity of EAS reactions. It is critical that prospective methods are easy for the end-users to work with to have a significant impact. This paper provides a fast, reliable, and easy-to-use computational tool to predict the site selectivity of EAS reactions more robustly than our previously reported NMR-based method. In a broader sense, the ability to block the more reactive site(s) in complex heteroaromatics may increase the use of halogens as protective groups in arene chemistry. It was recently reported that empirically calculated 1 H and 13 C chemical shifts can be applied to retrospectively account for the regiochemical outcome of ca 80% of 130 EAS reactions from the literature. For many of the remaining cases, the regioselectivity was rationalized by visual inspection of the HOMO orbital computed using DFT, bringing the success rate up to ca 95%. The chemical shifts were obtained immediately from the structure using ChemDraw. Conversely, DFT calculations are time-consuming for larger compounds and cannot be performed by non-expert users and there is a certain degree of subjectiveness in predicting the regioselectivity based on the HOMO orbitals. Finally, the sound judgment of chemists was required to evaluate the reactivity of unsubstituted and electron-deficient benzenes, and the method failed for 1,2,4-triazoles. In this paper the Kruszyk et al. study is extended to more than 525 reactions reported in the literature and we present a new semiempirical quantum mechanical (SQM)-based method that improves both the accuracy and precision of the predictions. The majority of the selected reactions were performed with N-bromosuccinimide (NBS) to avoid any substrate protonation that might occur when using for example bromine in acetic acid. This precaution may not have been necessary as the vast majority of cases where both conditions have been applied led to the same regiochemical outcome of the halogenation reaction. The computational work is based on the observation by Streitwieser and others that the rates of many EAS reactions correlate well with equilibrium values for protonation in solution. This observation implies that the site with the highest proton affinity (i.e. the protonated regioisomer with the lowest free energy) corresponds to the most probable site for EAS, as was demonstrated by Wang and Streitwieser for several polycyclic aromatic hydrocarbons. This approach is in line with the commonly accepted EAS reaction mechanism when considering the protonated species as a "surrogate" for the arenium ion (Scheme 1). The rationalization of the regioselectivities in the benzene series based on the relative stabilities of the possible sigma complexes is textbook knowledge. Jensen and co-workers recently showed that SQM methods can be used to accurately predict pK a values of ionizable groups. Gratifyingly, with only the SMILES string as input this approach pin-points the protonated regioisomer(s) with lowest energy and, hence, the most likely EAS reaction site(s).