I have been doing quite a lot of reading recently about the theory of atoms in molecules and think that I have semi-satisfactory answers to your questions, so I will share them here and perhaps others will add complementary answers as there actually seems to be quite an extensive literature on this subject, but there is still a lot of debate over these topics.
First, a quick note on how the charge density, $\rho$, can be used in the theory of atoms in molecules (AIM) to determine chemically intuitive properties. This is usually done by the definition of an index. For instance, the first index one might think to define is the electron number index, $N(\Omega)$
, where $\Omega$ is the space defined by the critical points surrounding a single nucleus. Essentially AIM claims that $\Omega$ and the very idea of an atom are interchangeable, so $N(\Omega)$ is conceptually identical to answering how many electrons are associated with a partiular atom. Hence, the AIM definition of the charge on an atom follows as,
$$
Q_{\Omega_i}=Z_{\Omega_i}-N(\Omega_i)
$$
where $Z_{\Omega_i}$ is the charge of the nucleus of atom $i$ in the basin $\Omega_i$.
Not surprisingly, $N(\Omega_i)$ is defined as[1],
$$
N(\Omega_i)=\int_{\Omega}\rho(\textbf{r})d\tau
$$
which is clearly just the charge density in a particular basin. This is interchangeable with a number density of electrons by division by the electric charge $e$.
These indices are always defined in terms of the charge density. The reason these are called indices, rather than something more exact, is because we are really computing a quantity which dimensionally satisfies the property we are interested in, and then interpreting it as the thing itself. The ambiguity comes from the idea of a basin, $\Omega$, which is mathematically well-defined but is not guaranteed to be physically exact for any reason.
How dependent are computed charges using the quantum theory of atoms in molecules on the used level of theory?
Ref. [2] is a review of the AIM theory with a specific emphasis on the electron localization functions (ELFs) to describe aromaticity. More on these ELFs in a moment. They present the calculation of the index $N(\Omega)$ for several diatomic molecules using DFT, HF, and CISD calculations. The total charge in each of the basins stays more or less the same for all of these methods. The value of $N(\Omega)$ generally only changes in the first or second decimal place for the calculations they present (so tenths to hundredths of an electron). The trend when Coulomb correlation is included is that $N(\Omega)$ decreases for what would traditionally be considered the more electronegative atom. That is to say that coulomb correlation makes the bonds more covalent and less ionic. For mostly ionic systems such as $\ce{LiF}$, the addition of coulomb correlation has a small effect because the electrons are already mostly localized in the basins [2]. All of this is to say that the actual charges computed based on $N(\Omega)$ seem to be relatively insensitive to coulomb correlation. This makes sense based on the fact that the primary features of the spatial distribution of $\rho$ will be dictated by the requirement of antisymmetrization of the wavefunction which is, of course, included in a HF description.
That covers the question of charges based on $N(\Omega)$, but this does not actually give us a picture of how electrons are shared between atoms. In other words, how localized and how delocalized are the electrons? In order to answer this question, ref. [3] provides definitons of both localization and delocalization indices, which are unique because the localization index describes localization of electrons to a single basin, while the delocalization index describes the sharing of electrons with either a specific basin, or with all other basin in a molecule.
While $N(\Omega)$ is not greatly affected by coulomb correlation, the localization index, $\lambda(A)$, and the delocalization index, $\delta(A,B)$ are. For instance, $\delta(A,B)$ which is a measure of the sharing of electrons between the atoms $A$ and $B$, changes by almost $1.0$ when going from the HF description of $\ce{N2}$ to the CISD description of $\ce{N2}$. In general, for homodiatomics, coulomb correlation has the effect of increasing the density in individual basins. And, consistent with the statement above about charge, the localization index for polar diatomics tends to increase the charge on the less electronegative element when a correlated wavefunction is used.
I would highly recommend reading ref. [3] as a place to gain some intuition about the behavior of these different indices.
Do calculated charges converge to a specific value with higher level of theory?
I have not been able to find any explicit discussion of this, but of course any method which approaches the exact wavefunction, and hence the exact density, must have the same limiting value for any index defined in terms of the density. This is really the power of the AIM model. It is completely quantum mechanical in the sense that the exact wavefunction will yield the exact basins. Whether or not these basins mean what AIM wants them to mean, however, is what people debate.
Another related point is how the various indices are affected by the particular basis set used. I have seen a couple papers which mention this, and note that the results really do not depend much on the particular basis set chosen[1]. The reason which is speculated for this behavior is that the partitioning used in AIM is in real space rather than in the Hilbert space to which the basis set belongs. Unsurprisingly, however, the indices are sensitive to the basis set when the density itself is very sensitive to the basis set. Open-shell molecules seem to be more complicated on this front [2].
Are charges from calculated electron densities comparable to charges from measured electron densities?
The primary purpose of ref. [1] is to compare the value of properties determined from AIM and from a density determined by x-ray diffraction of p-nitroaniline. This molecule was chosen, I believe, because it has a large dipole moment.
One complication in using the experimental density is that it is not as simple as just measuring the density and then having the whole scalar field in front of you. Rather, from what I understand, you measure the diffraction data and then use certain theoretical models to translate this diffraction data into a charge density. Thus, there is ambiguity in what charge density you actually measure. Ref. [1] uses the UMM and KRMM multipole formalisms which are discussed in this wikipedia article. The point of bringing this up is that the answer for the charges actually depends on which method you use to construct the experimental density.
Ref. [1] indicates that the KRMM multipole method which uses the AIM analysis of charges agrees much more closely with an AIM analysis of the charges from a theoretical density found from periodic-HF and periodic-DFT calculations (this is all solid phase). The p-HF and moreso the p-DFT charges agree very well with experimental charges which are calculated from KRMM using the AIM method. See Table 2 of [1].
The dipole moments also agree very well, but only after various refinements are made to both the partitioning of the experimental density as well as the use of the KRMM method. See Table 4 of Ref. [1].
One point which was mentioned in this book[4] (the relevant chapter is available on Google Books, ch. 3) is that the location of bond critical points, saddle points on the density between two atoms, can be very different for theoretical and experimental densities. Thus, properties derived from the basins can, in principle, be very different for the densities. It seems, however, that despite this problem of locating the bond critical point accurately, properties which are derived from the theoretical and experimental basins are not very different. This is because in the region of the saddle point, the density is usually quite flat, thus shifting the bond critical point does not have much effect on how much density is in each basin.
Are there more reliable methods to obtain partial charges?
Only briefly, I will point out that various population analyses are always possible to determine charges on atoms. The answer to this really depends on what you mean by reliable. In a sense, AIM is the only reliable method of determining charges on atoms in molecules because it is the only method I am aware of which bothers to define what an atom in a molecule actually is. Additionally, it is reliable in the sense that it exists in a completely self-contained mathematical framework. Population analyses, however, can be very dependent on the particular basis set used, and also are not invariant to unitary transformations of the orbitals. In this sense they are very unreliable.
Mulliken charges are a common form of finding charges from a population analysis, but these are known to very basis set dependent. A better scheme, I think, is the CHELPG scheme which fits the molecular electrostatic potential of a system and hence conserves the dipole moment when it assigns charges to atoms. This is a very nice property for a population analysis. Nonetheless, it will be sensitive to method and basis set depending on how much the electrostatic potential changes upon changing either of these.
I think that generally the AIM charges and charges from population analyses will not be too different, but one notable exception is given in ref. [4] on page 162 where it is noted that the AIM charge on boron in $\ce{(H3BNH3)2}$ is positive, $q(\Omega_{\ce{B}})=2.15$, while the Mulliken population analysis gives a negative charge of $-0.26$. So that's a very stark difference. I guess I would trust AIM, but I don't really know.
As a concluding point, everything I have discussed above calculates properties from indices based on $\rho$, but there is no reason to limit the definition of indices to the density. Rather, some authors have defined indices, for physically motivated reasons, based on the gradient of the density, $\nabla\rho$ and the Laplacian of the density, $\nabla^2\rho$.
An interesting application of one of these based on the gradient is for some alkanediols where there is no bond critical point between the two $\ce{O-H}$ groups despite the obvious presence of an internal hydrogen bond [5].
Ref. [5] is a really cool paper. Definitely give it a read.
[1] Volkov, A., Gatti, C., Abramov, Y., & Coppens, P. (2000). Evaluation of net atomic charges and atomic and molecular electrostatic moments through topological analysis of the experimental charge density. Acta Crystallographica Section A: Foundations of Crystallography, 56(3), 252-258.
[2] Poater, J., Duran, M., Sola, M., & Silvi, B. (2005). Theoretical evaluation of electron delocalization in aromatic molecules by means of atoms in molecules (AIM) and electron localization function (ELF) topological approaches. Chemical reviews, 105(10), 3911-3947.
[3] Fradera, X., Austen, M. A., & Bader, R. F. (1999). The Lewis model and beyond. The Journal of Physical Chemistry A, 103(2), 304-314.
[4] Popelier, P. L. A., Aicken, F. M., & O’Brien, S. E. (2000). Atoms in molecules. Chemical Modelling: Applications and Theory, 1, 143À198.
[5] Lane, J. R., Contreras-García, J., Piquemal, J. P., Miller, B. J., & Kjaergaard, H. G. (2013). Are bond critical points really critical for hydrogen bonding?. Journal of chemical theory and computation, 9(8), 3263-3266.
