Getting the forces on Silicon with a 2x1x1 kgrid gets a nonzero result (irrespective of symmetries settings). This is consistent with finite differences. OTOH abinit gives zero by default. Setting nsym to 1 in abinit (ie deactivating the symmetries) does give a nonzero result. Basically what's happening is that abinit symmetrizes everything wrt the model symmetries, and we symmetrize wrt the basis symmetries. I like what we do better, because the symmetries are then a pure optimization (unlike abinit, where switching off symmetries gives a different result, even for the ground state energy). Not sure what other codes do. We should document it, though.