An explicit application of isodesmic reaction (a proton exchange between the studied and similar in structure reference molecule), where the free energy change of the protonation reaction in water was obtained using the free energies in solution from a single continuum model, was used to predict stepwise protonation constants of nitrilotriacetic acid.Calculations were performed at the RB3LYP/6-311+G(d,p) level of theory in conjunction with PCM-UA0 solvation model. Five reference molecules were investigated. It has been established that one must pay a special attention to structural similarities between the studied and reference molecules and selection of a protonated form of the reference molecule. The protonation reactions in which the studied and reference molecule are involved in must be (if possible) of the same order; e.g. first (or generally nth) protonation reaction of the reference molecule must be used to compute the first (or nth) protonation constant of studied molecule. The lowest energy conformer must always be used. The first, second, third and fourth computed protonation constants differed, on average, from experimental values by 3.3, 0.8, 0.2 and 0.2 log units, respectively. It appears that the charge on the reference molecule has more decisive influence on accuracy of computed protonation constants than its structural differences when compared with the studied molecule. Results reported can be used as guide in constructing isodesmic reactions useful for the theoretical prediction of protonation constants by use of methodology described in this work.-2 -