Without any constraint the symmetry of the Cr 3d orbital occupations will break and Cr atoms equivalent by symmetry will not show exactly the same occupation pattern on convergence. The numerical noise introduced by the employed grids in CP2K will always trigger such a break in symmetry, even for insanely high cutoffs and already small fluctuations in the 3d occupation patterns of topological equivalent Cr atoms can easily cause energy differences up to a milliHartree on convergence as you observe. CP2K has currently no tools implemented to impose the required symmetry strictly. I tried to keep at least the 3d occupations of the Cr atoms spherical by using an AMF setting which results in the following energies for the system configuration you provided (GPW, 2000 Ry cutoff, EPS_DEFAULT 1.0E-16):
# System Energy [a.u.]
The total 3d occupations of symmetry equivalent Cr atoms, however, can still vary slightly, which causes the microHartree fluctuations in the total energy. I am afraid, it will be difficult to get it much better.
I enforced a smearing of the Cr 3d orbital occupations by adding the following section in the &DFT_PLUS_U of each Cr kind:
ORBITALS -2 -1 +0 +1 +2
Note, with that, you do not converge to the ground state. This would require the implementation of a symmetry constraint which keeps the same 3d occupation pattern for all Cr atoms equivalent by symmetry.
You can assign each group of Cr ions equivalent by symmetry to a specific atomic kind and then enforce the desired d orbital occupations for that atomic kind and hope that the initial symmetry won’t break during the wavefunction optimization. Currently, CP2K provides no measure to keep an equal number of d electrons (smeared within the selected subset of d orbitals) per symmetry equivalent Cr ion within such a group during the SCF procedure. You will need to hack the code to achieve that additional symmetry constraint. You can contact me directly when you need help for that.
It is always good to be optimistic. I agree with you that it might be well possible to achieve the desired energy resolution for relative energies with and without perturbation in static calculations (i.e. without an ongoing cell size or shape change or significant atomic movements).