Quantum chemical calculations were carried out for deprotonated (P−) and protonated purine (PH+) and for adducts with one alkali metal cation (P−M+ and PM+, where M+ is Li+ or Na+) in the gas phase {B3LYP/6-311+G(d,p)}, a model of perfectly apolar environment, and for selected structures in aqueous solution {PCM(water)//B3LYP/6-311+G(d,p)}, a reference polar medium for biological studies. All potential isomers of purine derivatives were considered, the favored structures indicated, and the preferred sites for protonation/deprotonation and cationization reactions determined. Proton and metal cation basicities of purine in the gas phase were discussed and compared with those of imidazole and pyrimidine. Bond-length alternations in the P, PH+, P−M+, and PM+ forms were quantitatively measured using the harmonic oscillator model of electron delocalization (HOMED) indices and compared with those for P. Variations of the HOMED values when proceeding from the purine structural building blocks, pyrimidine and imidazole, to the bicyclic purine system were also examined. Generally, the isolated NH isomers exhibit a strongly delocalized π-system (HOMED > 0.8). Deprotonation slightly increases the HOMED values, whereas protonation and cationization change the HOMED indices in different way. For bidentate M+-adducts, the HOMED values are larger than 0.9 like for the largely delocalized P−. The HOMED values correlate well in a comprehensive relationship with the relative Gibbs energies (ΔG) calculated for individual isomers whatever the purine form is, neutral, protonated, or cationized. When PCM-DFT model was utilized for P−, PH+, PM+, and P−M+ (M+ = Li+) both electron delocalization and relative stability are different from those for the molecules in vacuo. The solvation effects cause a slight increase in HOMEDs, whereas the ΔEs decrease, but in different ways. Hence, contribution of particular isomers in the isomeric mixtures of PH+, PM+, and P−M+ also varies.