The complete 1-loop computation of the processes $(\tchi^0_i \tchi^0_j \leftrightarrow����, ~ gg)$, for any pair of the four MSSM neutralinos, has been performed for an arbitrary c.m energy. As a first application suitable for Dark Matter (DM) searches, the neutralino-neutralino annihilation is studied at relative velocities $v/c\simeq 10^{-3}$ describing the present DM distribution in the galactic halo, and $v/c\simeq 0.5$ determining neutralino relic density contributions. Numerical results are presented for 31 MSSM benchmark models indicating considerable sensitivity to the input parameters. Numerical codes are released, which may be used for the computation of the annihilation of any two neutralinos to two photons or two gluons at the aforementioned $v/c$ values. In the near future, we intend to complement them with codes for the inverse process $����\to \tchi^0_i \tchi^0_j$, observable at the future high energy Linear Colliders.

14 pages and 1 figure; version to appear in Phys. Rev.D; Title slightly changed