Check2xsf Visualisation: Notes
A couple of points which it may be worth considering.
If ultrasoft pseudopotentials are used, the wavefunctions are not normalised, and the total charge density is not the sum of the individual wavefunctions squared. This is a feature of ultrasofts, not a bug...
If density mixing is used, the total charge density is the sum of the individual wavefunctions squared only when the wavefunctions are fully converged.