Hi,
two things you could/might have to do:
1. grid
To catch all of the density you probably need to use a bigger box which reaches further out in space - depending on the basis set.
(a bit off topic):
To find out which atomic orbitals contribute to the density (by checking the occupied molecular orbitals and printing the contribution of each atomic orbital) try:
tm2molden mostat <list-of-MOs> above <threshold>
So if you have 40 occupied orbitals (eiger quickly tells you how many you have) and you want to see all AOs which contribute to each MO more than, say, 0.01% just call on the command line:
tm2molden mostat 1-40 above 0.01
Enter a few times <enter> and then you will get a file named mostat.out which looks like that:
MO Index (from eiger): 3
Sym= 3a
Energy [Hartree]= -.99177169241775E+01
Energy [eV]= -.26987482366088E+03
Spin= Alpha
Occup= 2.000000
contribution of AO 2 c 1 s is: 24.939 %
contribution of AO 6 c 1 s is: 24.939 %
contribution of AO 3 c 1 s is: 24.939 %
contribution of AO 5 c 1 s is: 24.939 %
contribution of AO 2 c 2 s is: 0.052 %
and so on. That is useful to see how many diffuse AOs contribute to the occupied orbitals and thus to the density.
2. numerical accuracy
Another source of numerical inaccuracy is the format of either plt or cube files. They are by default single-precision data. As the density is written to a grid and then summed up, just from pure numerics the sum will never exactly be the number of electrons (well, never is not correct, it could coincidentally be just an integer number...).
Try to add fmt=dtx to the $pointval input line instead. This will print double precision data as ASCII text (dtx = double precision text file). Looks like that:
#origin 0.000000 0.000000 0.000000
#vector1 1.000000 0.000000 0.000000
#vector2 0.000000 1.000000 0.000000
#vector3 0.000000 0.000000 1.000000
#grid1 start -5.767470 delta 0.202135 points 66
#grid2 start -5.767470 delta 0.202829 points 59
#grid3 start -5.020450 delta 0.204532 points 67
#title for this grid 101
#density
#plotdata
# cartesian coordinates x,y,z and f(x,y,z)
-5.7674700000000E+00 -5.7674700000000E+00 -5.0204500000000E+00 3.9793580979130E-10
-5.5653350769231E+00 -5.7674700000000E+00 -5.0204500000000E+00 5.8733532389287E-10
-5.3632001538462E+00 -5.7674700000000E+00 -5.0204500000000E+00 8.5021125663768E-10
-5.1610652307692E+00 -5.7674700000000E+00 -5.0204500000000E+00 1.2071193363237E-09
...
The last column then needs to be summed up. Whether or not that will increase accuracy - it's at least worth a try.