  We describe low-lying collective excitations of atomic nuclei with the
multi-reference covariant density functional theory, and combine them with
coupled-channels calculations for heavy-ion fusion reactions at energies around
the Coulomb barrier. To this end, we use the calculated transition strengths
among several collective states as inputs to the coupled-channels calculations.
This approach provides a natural way to describe anharmonic multi-phonon
excitations as well as a deviation of rotational excitations from a simple
rigid rotor. We apply this method to subbarrier fusion reactions of
$^{58}$Ni+$^{58}$Ni, $^{58}$Ni+$^{60}$Ni and $^{40}$Ca+$^{58}$Ni systems. We
find that the effect of anharmonicity tends to smear the fusion barrier
distributions, better reproducing the experimental data compared to the
calculations in the harmonic oscillator limit.
