It has been demonstrated for the first time that prediction of several consecutive protonation constants for the highly and negatively charged molecules, such as nitrilotripropanoic acid (NTPA), is possible with acceptable accuracy when isodesmic reaction (IRn) methodology, instead of commonly employed thermodynamic cycle (TC), is employed. Four stepwise protonation constants of NTPA were computed (RB3LYP/6-311+G(d,p) level of theory employing PCM/UA0 solvation model) to within 1 log unit of experimental data with an average error in the protonation constant of about 0.5 log unit. This good agreement was achieved for minimum energy structures of NTPA (studied ligand) and iminodiacetic acid (reference molecule). Results obtained strongly support the view that full conformational analysis should be seen as prerequisite for computing protonation/dissociation constants from IRn and possibly also from TC. Methodology proposed here broadens up, in our opinion, a scope of studying protonation constants computationally and opens up a new field of applications for poly charged ligands. TC did not work here at all as proton on N-atom was not preserved in gas-optimized structures; this proton always protonated available COO-group instead.