This is a really annoying problem! I believe it happens because of finite numeric thresholds and most often we observe it when there are degenerate orbitals, like in linear molecules or in the noninteracting fragments that you described.
Because any linear combination of degenerate orbitals is valid, the SCF procedure can mix them, unless there are special controls. I also think it happens more often with large bases with diffuse orbitals.
There are several things to try:
 Use new integral library by setting ‘use_libqints=true’, it shows a more robust behavior.
 Do a fake DFT calculation first (using the grid helps to resolve symmetry
of degenerate irreps):
$rem
…
exchange gen
…
$end
$xc_functional
X HF 1.0
X Slater 0.000001
$end
and then read the orbitals into the second job using ‘scf_guess = read’

Sometimes playing with guess and SCF algorithms can help.

If desperate, try to use a smaller basis.