The optical counterpart of the gravitational wave event GW170817, known as kilonova, has provided strong evidence that binary neutron star mergers are favourable sites to host the r-process nucleosynthesis. Kilonova is a quasi-thermal electromagnetic emission powered by the radioactive decay of heavy neutron-rich nuclei produced by the r-process. Considering the variety of elements contributing to kilonova ejecta, essential information about its composition can be achieved through spectral characterisation, radiative transfer simulations, and opacities. The latter represents one of the most challenging aspects of the modelling, as it relies on accurate atomic structure calculations of energy levels and transitions. Since light r-process elements are major opacity contributors in early (< 2 days) scenario, this work focuses on atomic calculations for Zr I-IV. Energy levels and bound-bound transitions are determined using the GRASP2018 code, assuming two different multi-reference sets for each ionisation stage: one including, and one excluding core-core and core-valence correlations. Results demonstrate that the inclusion of f shell and core correlations impacts on both energy levels and transitions. A systematic assessment of the accuracy is performed through detailed comparisons with the NIST ASD. Finally, these Zr data are integrated on the open access MARTINI platform.