The aim of this research was to study the stability of the structure of the ZrO2 doped with trivalent oxide Zr1-xMxO2-δ (M = La3+, Nd3+, Sm3+, Eu3+, Gd3+, Y3+, Er3+, Yb3+ and Lu3+ through atomistic modelling and bond valence sum method. Short range potential used in this study was Buckinghams’ potential. Result of geometry optimization at constant pressure shown both cell parameters of ZrO2 was in good agreement with experimental results because of the difference was only 0.11%. Increasing the concentration and the size of substituting dopant of ZrO2 makes the lattice energy of the doped structure was more positive so that the stability of the doped ZrO2 structure decreases. The decrease in the stability of ZrO2 doped with Y3+, Er3+, Yb3+ and Lu3+was smaller than ZrO2 doped with La3+, Nd3+, Sm3+, Eu3+ and Gd3+. BVS results shown that the structure of ZrO2 doped with La3+was not appropriate because it has different value of BVS was more than 0.1