Kinematic error is an important indicator for evaluating the performance of harmonic drive gears. Previous methods for calculating the kinematic error only considered a single meshing tooth, the effect of other meshing teeth was compensated by the error homogenization coefficient. Because of the complex contact state of harmonic drive gears, the results of these methods are always inconsistent with experimental results. This paper proposes a mathematical model considering all meshing teeth to predict the pure kinematic error of harmonic drive gears. First, the contact points on each meshing tooth are determined, and the contact states of the flexspline and circular spline teeth are investigated. Then, the contact forces on the meshing teeth are obtained, and the pure kinematic error is calculated based on the relationship between the torque caused by the contact force on the left and the right tooth profiles. Finally, the mathematical model is validated by simulation and experimental methods. The study in this paper can not only be used to accurately predict the pure kinematic error but also is valuable in the selective assembly of higher precision harmonic drive gears.