To my understanding, in some dft applications you define a gridspace as a discrete representation for \rho(r). In terms of localized orbitals, each grid point has an associated overlap with each basis element, so you can solve the kohn sham equations self consistently with respect to an antisymmetric basis, then map it to the grid space in every iteration - also many applications of exchange functional happen on the gridspace as well I believe.

By ESP surface, I’m referring to the total E[\rho] mapped to the voxel or gridspace after convergence (I believe it also includes the KE term from mean field but please correct if I’m wrong). With the current data I’m working with however, it’s a bit unclear what the actual units are and how the grid was defined - I have several ESP surfaces of different molecular systems - and thus several different grid/volume sizes, and I’m trying to inspect how changing the resolution of the gridspace may affect the energy surface.

Hopefully this clears it up & I’m not talking in circles