  We develop a rigorous and highly accurate technique for 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 which 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 two-fold 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 on the order of $10^{-5}$ or lower. We
observe exponential decay of the error with the distance between vortices,
studied in the range from one to nine plaquettes. We also study its correlation
with the involved energy gaps and provide preliminary analysis of the relevant
adiabaticity conditions. The work allows to investigate the Berry phase in
other lattice models including the Yao-Kivelson model and particularly the
square-octagon model. It also opens the possibility of studying the Berry phase
under non-adiabatic and other effects which may constitute important sources of
errors in topological quantum computation.
