The electron density of a molecule is very far from an equally distributed charge density. Therefore, the common density sampling on a regular grid for visualisation purposes is not useful for cases where numerical properties of the density distribution are sought after. In that case, integration grids as Becke-Lebedev grids are required. Based on the atomic coordinates, they give a set of points in space with associated weights to account for the unequal volume represented by each sampling point in the space-filling partition of the electron density. Now the problem is to evaluate the electron density at the points in question rather than on a regular rectangular grid.
In the context of Gaussian calculations, here is how. The feature is implemented in the tool cubegen albeit it lacks clear documentation.
First, create a text file (e.g. gridpoints) with the points at which you’d like to sample the electron density. Each row should contain the space-separated coordinates of a single point in angstroms. Now run the calculation and create a formatted checkpoint file (say run.fchk) from the checkpoint file of your calculation with formchk. With the command
cat pts | cubegen 1 Density=CC run.fchk test.cube -5
you get a new file, test.cube, which contains your input points and as a fourth column the density in atomic units. While the extension hints otherwise, this is not a valid cube file, but is a machine readable representation of the electron density on the grid you have specified. Note that typically you still need the weights of the individual grid points. The code that generates the grid point coordinates for you should also give you the weights.
Using a non-rectangular grid is also faster for integration purposes and has a lower memory footprint since the number of grid points for a given quality of the result is much smaller. Note that rectangular grids are much better suited at visualisation purposes than the integration grids. Always try to use tools for the purpose they’re designed for…