Here are a few notes that I completely forgot when answered questions:
About the role of weak correlation: it is extremely hard to assess it, because about 95% of the wave function is the part that the model space contain (see its definition below). Only very small part of the wave function contains other determinants that are not included in the model space directly - they are wrapped in through the effective interactions. These determinants are the cause of the weak correlation. The overall change in the gap is the differential effect - the weak correlation of the triplet and the singlet are slightly different. I am not aware of a rigorous study that would investigate that, but some simplistic models predict that this could be due to a tiny fraction of charge transfer between ligands and metal centers. I could be wrong here - please feel free to correct me.
About spin contamination. Basically, the extraction of the triplet density matrix works in AO - spin contamination does not prevent the algorithm to work. The only issue is that one needs to assign multiplicity to the state to have proper dimensions of the Hamiltonian matrix and to evaluate the correct Clebsh-Gordan coefficients. When the spin contamination is not large, correct <S^2> is easy to guess. When spin contamination is very large, it is not possible to tell much, because it is not possible to say what would be the multiplicity. The algorithm does not work for broken spin case, for example.
About broken-spin DFT and orbitals. Let me first say what we do for the orbital localization. We compute UHF for the high-spin state, identify open-shell orbitals, and proceed to CCSD and EOM-SF-CCSD. The open-shell alpha orbitals are found through SVD of alpha-beta occupied-virtual overlap matrix. There are very few singular values close to one (the corresponding singular vectors are the alpha and beta open-shell orbitals), while the remaining singular values are close to zero. The latter may not be zero is spin contamination is presented - this directly translates to <S^2> expression for UHF. Open-shell orbitals lead to a well-defines localization problem if they sit on different centers, as is the case of the copper systems that I considered. These orbitals contain some contributions of ligands. The wave functions are very compact in these orbitals. Canonical orbitals are very delocalized and messy; it is hard to make physical sense from them. The EOM amplitudes in these orbitals contains hundreds of important weights - the wave function is not compact.
Natural orbitals seem to be a bit more compact than the open-shell orbitals. This could be due to a partial correlated screening that the open-shell electron experiences around adjacent ligands.
Broken-spin DFT is not the most reliable method. Sometimes it works, sometimes it does not. In EOM-CC we, at least, know when EOM-CC fail - we always monitor which contributions come into the wave function and based on that we judge the quality of the calculation. Broken-spin DFT lacks this ability. The character of the orbital depends on how the orbital was selected for localization. For KS-DFT I see that the KS orbitals have the same problem as the HF canonical orbitals - they are too delocalized and messy. I would recommend the overlap SVD procedure for broken-spin DFT first to see whether the problem remains.