H()

After the above modifications, the only remaining call of cI() on a wave-function object is in the application of the local potential in H(), in the form cIdag(Diagprod(gbl_Vdual+Vsc,cI(C))). Rewrite this contribution to H() using an explicit loop over the columns of C with a code fragment along the lines of

for col=1:size(C,2)
 ... cIdag((gbl_Vdual+Vsc).*cI(C(:,col))) ...
end

Note: You may have used a different internal variable name for Vsc, the sum of the nuclear, Hartree and exchange-correlation potentials.

Debugging: Rerun ``Hatoms.m'' to convince yourself that the output is unchanged.