Concerning your specific question, I think a direct application of Kostant's multiplicity formula will show that the multiplicity of the zero weight in the finite dimensional module with highest weight $r\rho$ for $r>0$ is always given by applying the partition function to $r\rho$. (I'm not sure what symbol you are using, but at one point it reverts to a familiar $P$.)
In more detail, $\rho$ has two basic characterizations: as the sum of fundamental dominant weights (showing that it is dominant) and as the half-sum of all positive roots (showing at least that $2\rho$ lies in the root lattice, though $\rho$ itself might or might not). From representation theory it follows that the value of the partition function at $r\rho$ is the dimension of the zero weight space of the corresponding infinite dimensional Verma module.
Kostant's basic insight was that the finite dimensional simple quotient of any Verma module with dominant highest weight has formal character given by an alternating sum of naturally related Verma module characters (which later got encoded in the BGG resolution). So the underlying representation theory problem is to show that nothing in the alternating sum in this special case will reduce the size of the zero weight space. That's a direct consequence of Kostant's formula, since the other weights involved (with Weyl group action shifted by $\rho$) aren't dominant in the positive root cone when you are computing just the multiplicity of the zero weight.
More generally, Kostant explained in the introduction to his 1959 paper here how his intuition developed toward the rigorous formula. Later work by Bernstein-Gelfand-Gelfand streamlined the derivation of this formula (and Cartier pointed out earlier the logical equivalence of the formula with Weyl's), but it may be useful to take a look back at the original paper. [Note that the symbol $g$ is used there for what later became standardized as $\rho$.]
Concerning asymptotics of the partition function, I'm not sure how much has been written down. But for example an appendix in Verma's 1966 Yale thesis makes some observations along this line; as I recall he developed a polynomial bound. Much later Lusztig (and then R.K. Brylinski) studied a $q$-version of the partition function as well. Such work is done in a purely theoretical mode, of course, but there may be other computational approaches of interest for your question.
P.S. Sorry for the imprecise wording in my earlier version. To get a nonzero value for the partition function you need an element in the positive root cone
($\mathbb{Z}^+$-linear combination of simple roots). But any $w \neq 1$ in the Weyl group changes the sign of some simple root in a multiple of $\rho$, after which subtracting $\rho$ makes it even worse. (Here and in the asymptotics, it's best to express elements of the root lattice uniquely as combinations of simple roots; then the rank of the root system plays a natural role in any asymptotic estimate.)
