An explicit application of isodesmic reaction (a proton exchange between the Studied and Structurally similar 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 the PCM-UA0 solvation model. Five reference molecules were investigated. It has been established that one must pay 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., the first (or generally nth) protonation reaction of the reference molecule must be used to compute the first (or nth) protonation constant of the Studied molecule. The lowest energy conformer must always be used. The first, second, third, and fourth computed protonation constants differed, oil 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 the accuracy of computed protonation constants than its structural differences when compared with the studied Molecule. Results reported call be used as a guide ill Constructing isodesmic reactions useful for the theoretical prediction of protonation constants by use of methodology described in this work.