intervoxel said:
1) Yes, I'm using Griffiths formulae, which uses the SI system.
2) I followed the recommendation contained in the paper by Tarantola: http://www.cmmacs.ernet.in/inverse/Diverse/DensityVersusVolume/DensityVersusVolume.pdf (equation 9).
[STRIKE]Ok, that looks bizarre and wrong, at least for applications in Q.M. Everything in section 3 is complete nonsense if one is using the standard definition of a probability density in Q.M., which is just the square modulus of the wavefunction multiplied by the volume element.[/STRIKE] You must multiply by the proper volume element ... r^{2}sin\theta drd\theta d\phi ... that is part of the definition of the probability density in Q.M.
EDIT: Ok, I wrote the crossed out stuff above too quickly .. after looking at the paper more carefully, everything the author wrote is correct ... his phrasing is just a little weird, and section 3 really threw me for a loop when I was skimming through.
3) Nope! I modified the program but the problem still remains (I divided the integrated value by the number of elements).
step 3 shouldn't be necessary *IF* you are using the correctly scaled volume element for your sum (including units).
EDIT: Think about it this way, what is the volume of a region of space bounded by equal sized steps along the coordinates r, theta and phi? The answer depends on where you are relative to the origin right? Regions close to the origin will have smaller volumes than those that are farther away, and you must account for this somehow. You can either do it by using a probability density, which includes the volume element in the definition (this is what the author is talking about in his paper). Or you can just integrate \psi^{*}\psi, which is what the author calls a "volumetric probability", but then you need to account for the volume element explicitly when doing the integration. Both treatments end up the same for your case:
P\left(r,\theta,\phi\right)=\int^{2\pi}_{0}d\phi\int^{\pi}_{0}sin\theta d\theta\int^{\infty}_{0}r^{2}dr\psi^{*}\left(r,\theta,\phi\right)\psi\left(r,\theta,\phi\right) \approx \sum_{i=1}^{M}r^{2}_{i}\Delta r\sum_{j=1}^{N}sin\theta_{j}\Delta\theta\sum_{k=1}^{Q}\Delta\phi\left[\psi^{*}\left(r_{i},\theta_{j},\phi_{k}\right)\psi\left(r_{i},\theta_{j},\phi_{k}\right)\right]
where all points along a given coordinate are equally spaced, and N\Delta\theta=\pi, Q\Delta\phi=2\pi, M\Delta r=R
where R is some large arbitrary value, like your choice of 30 Angstroms.