The prediction of conformations and protonation sites for the macrocyclic ligands H₂DO2A (1,4,7,10-tetraazacyclododecane-1,7-diacetic acid) and H₂ODO2A (1-oxa-4,7,10-triazacyclododecane-4,10-diacetic acid) has been performed employing the simulated annealing (SA) method and density functional theory (DFT) calculations using the B3LYP/6-31G* method in a vacuum and aqueous solution. These SA method/DFT calculations reveal that, in contrast to the H₂ODO2A ligand system, the H₂DO2A ligand system is (i) pre-organized for trivalent lanthanide (Ln) and other metal ion complexation, (ii) structurally more symmetrical and slightly more compact in aqueous solution (i.e. more and/or shorter intra-molecular hydrogen bonds), and (iii) with a greater degree of partial positive charge accumulation on the hydrogen atoms bonded to macrocyclic ring nitrogen atoms when protonated. The H₂ODO2A ligand system is not pre-organized. These observations are in accord with the experimental findings that the LnODO2A⁺ complexes are less thermodynamically stable and kinetically more labile as compared to those of the corresponding LnDO2A⁺ complexes. The results on the prediction of the ligand protonation sites are consistent with those experimentally obtained via NMR spectroscopy. The calculations of the first and second protonation constants are, however, not as accurate as compared to those experimentally determined using either the thermodynamic cycle (TC) or the isodesmic reaction (IRn) methodology, although the latter gave relatively better results. The lowest energy structures of the LnL⁺ and ZnL (Ln = Eu, Y; L = DO2A, ODO2A) complexes are also calculated using the same method. The Gibb's free energies (ΔG_aq) for a number of ligand and/or metal ion exchange reactions such as LnDO2A⁺ + H_nODO2A⁽²⁻ⁿ⁾⁻ LnODO2A⁺ + H_nDO2A⁽²⁻ⁿ⁾⁻ (Ln = Eu, Y; n = 0, 1, 2), LnDO2A⁺ + Ln′ODO2A⁺ LnODO2A⁺ + Ln′DO2A⁺ (Ln = Eu, Ln′ = Y), LnDO2A⁺ + Ln′³⁺ Ln′ODO2A⁺ + Ln³⁺ (Ln = Eu, Ln′ = Y) and LnDO2A⁺ + ZnODO2A LnODO2A⁺ + ZnDO2A (Ln = Eu, Y) have been calculated in aqueous phase and the reaction directions in some cases could be predicted to be consistent with experimental or expected results. The errors between the calculated and experimental Gibb's free energy data are in the range ΔG_aq,calc − ΔG_aq,exp = −1.89 to +7.00 kcal mol⁻¹ in seven selected cases involving LnDO2A⁺, LnODO2A⁺ (Ln = Eu, Y), ZnDO2A and ZnODO2A complexes. The predicted reaction directions with the small core effective core potential (ECP) data are not necessarily better than those using large core ECP. However, the former takes much longer computer time to obtain the energy data.