Adaptive mesh refinement (AMR) has been shown to allow solving partial differential equations to significantly higher accuracy at reduced numerical cost. This paper presents a state-of-the-art AMR algorithm applied to the multigroup neutron diffusion equation for reactor applications. In order to follow the physics closely, energy group-dependent meshes are employed. We present a novel algorithm for assembling the terms coupling shape functions from different meshes and show how it can be made efficient by deriving all meshes from a common coarse mesh by hierarchic refinement. Our methods are formulated using conforming finite elements of any order, for any number of energy groups. The spatial error distribution is assessed with a generalization of an error estimator originally derived for the Poisson equation.Our implementation of this algorithm is based on the widely used Open Source adaptive finite element library deal.II and is made available as part of this library's extensively documented tutorial. We illustrate our methods with results for 2-D and 3-D reactor simulations using 2 and 7 energy groups, and using conforming finite elements of polynomial degree up to 6.