Occupancies

The various energy and gradient functions require access to the filling filling factor $f$. Rather than passing $f$ as an argument down through all of the various functions, we will allow ourselves one final global variable, ``gbl_f'', a blemish which a C++ implementation could easily eliminate. Add the following lines to the very top of ``dft.m'',

%# Set the orbital occupancies
global gbl_f;
gbl_f=2; %# The usual case of a constant filling of two electrons per orbital

Your original Schrödinger solver implicitly assumed state occupancies of $f=1$. The above energy and gradient formulas for density-functional theory should work just as well for general values of $f$, even if the $\phi$ and $\epsilon_{xc}$ terms are set to zero. Thus, before adding any of the new terms, add the factors of $f$ to the appropriate places in getE() and getgrad().
Hints: Don't forget to declare gbl_f as a global variable in each function which needs ``f'' and to include it in both the density and the kinetic energy parts.