We develop a rigorous and highly accurate technique for the calculation of the Berry phase in systems with a quadratic Hamiltonian within the context of the Kitaev honeycomb lattice model. The method is based on the recently found solution of the model that uses the Jordan–Wigner-type fermionization in an exact effective spin-hardcore boson representation. We specifically simulate the braiding of two non-Abelian vortices (anyons) in a four-vortex system characterized by a twofold degenerate ground state. The result of the braiding is the non-Abelian Berry matrix, which is in excellent agreement with the predictions of the effective field theory. The most precise results of our simulation are characterized by an error of the order of 10−5 or lower.We observe exponential decay of the error with the distance between vortices, studied in the range of one to nine plaquettes. We also study its correlation with the involved energy gaps and provide a preliminary analysis of the relevant adiabaticity conditions. The work allows one to investigate the Berry phase in other lattice models including the Yao–Kivelson model and particularly the square–octagon model. It also opens up the possibility of studying the Berry phase under nonadiabatic and other effects that may constitute important sources of errors in topological quantum computation.