Given a complete separable metric space $(X,d)$. My question is if it is true that for any finite positive Borel measure $\mu\in\mathcal{M}^+_b(X)$ there exists a sequence of finite atomic measures $a_n\in\mathcal M^+_b(X)$ such that $$ a_n\rightharpoonup \mu, $$ i.e. $a_n$ converges narrowly to $\mu$.
The most obvious guess is that given $x_n$ a dense sequence in $X$ and $U^i_n$ are disjoint open sets such that $$ \bigcup_i U^i_n=B(x_0,n) $$ we define $$ a_n=\sum_{i=1}^{g(n)}\delta_{x_i}\mu(U^i_n), $$ where $g$ should have the property $$ \lim_n \frac{g(n)}{\mu(B(x_0,n))}=C. $$
Then we should have that for any continuous bounded function we have $$ \int_X f\ da_n=\sum_{i=1}^{g(n)} f(x_i)\mu(U^i_n)\rightarrow\int_Xf\ d\mu. $$ It should happens because $$ \sum_{i=1}^{g(n)} \min_{U_n^i}(f)\chi_{U^i_n}\leq f\leq\sum_{i=1}^{g(n)} \max_{U_n^i}(f)\chi_{U^i_n} $$ but there is no reason to infer that the choice of $U_n^i$ is good enough to refine the step functions from above and below well enough to obtain the integrals. Because we would need the choice to be independent of the function $f$.
What am I missing? Is it possible to begin with? Where can I find this result in the literature?
That's not a general answer, but it works in particular cases.
If $\mu$ is a probability measure and $X$ is compact (so that the set of continuous bounded functions is separable), there is an elegant way to show the result. Take $f_n$ a dense sequences of functions, and let $Y_1,\ldots,$ be iid in $X$ according to $\mu$. Then, for each $n$, almost surely, $\frac{1}{p}\sum_{k=1}^p{f_n(Y_k)} \rightarrow \int_X{f_nd\mu}$.
So, for almost every infinite sequence $x_n$ of elements of $X$, for each $p$, $\int_X{f_pd\mu_n} \rightarrow \int_X{f_pd\mu}$, where $\mu_n=\frac{1}{n}\sum_{k=1}^n{\delta_{x_k}}$. This implies that $\mu_n$ converges narrowly to $\mu$.
I think that you can still manage something when the measure has compact support (by restricting to the compact case). So in the case of Polish spaces, where finite measures are supported "up to $\epsilon$" on compact subsets, it may be possible to manage something as well.