Hello Yike,
Basically, the Harris-Foulkes functional can be seen as first order approximation to the Kohn-Sham functional:
Let's say you know the ground-state density up to a certain delta-rho , so that n^out = n^in + delta-rho.
Now you can replace n^out in the KS-functional, expand the expression, and truncate the resulting second order terms in delta-rho.
This results in the appearance of the Coulomb and XC "correction" terms which you are seeing,
and - more importantly - in an energy functional expression which is only dependent on the input density.
A good starting point is the book "Interatomic forces in condensed matter" by Mike Finnis.
Kind regards,
Fabian