The ``gmx density -dens number’’ gives a number density of around 33 water molecules per \text{nm}^3 (along any direction of simulation box) which perfectly agrees with what one gets if calculates the number of water molecule in a \text{nm}^3 just by considering the \text{MW}_{\text{water}} = 18 \frac{\text{g}}{\text{mol}} and \rho_{\text{water}} \approx 1 \frac{\text{g}}{\text{cm}^3} and Avogadro number.
However, I am not sure if the ``gmx density -dens number’’ works fine for other molecules. I just did the same as water above for another small molecule, I get a correct mass density for that molecule, however, the number density (\frac{\text{number of molecule}}{\text{nm}^3}) is much higher than what I get when I use the \text{MW} and \rho of the molecule and Avogadro number to calculate the number of the molecule that accommodate in one \text{nm}^3 of volume.
I wonder if that is possible, or I am doing something wrong in my calculations?
I guess GROMACS might be calculating number density by just dividing number of molecules by volume of the simulation box.
Few things you can cross check are:
Check whether you get the same answer by getting volume from gmx energy and then calculating number density by dividing total number of molecules by volume. If there is significant difference then one possible error could be you have not saved enough trajectories.
Check whether you are calculating density in properly equilibrated trajectory.
Check the phase of other molecules you are testing, water is in liquid phase, but for your molecules, for example, if they are in solid/gel phase, they may take time to equilibrate, or they may not equilibrate at all, and just vibrate around initial positions given in initial gro file. Here you can use annealing etc. for speeding up the dynamics to achieve equilibration.
There could be other reasons, but I would suggest first check above suggestions.