Here are a few clues that might help.

Note that the $L^3$ conjugated exponent is $3'=3/2$. For given $R>0$ let $B_R=B_R(0)$ be the ball of radius $R$ centered at the origin. Using spherical coordinates with $dx\sim r^2 dr$ in dimension $3$ you can compute explicitly
$$
\left|\frac{1}{|x|}\right|_{L^{3/2}(B_R)}=CR
$$
for some universal $C>0$.

**Claim 1:** if $\mathcal{A}=\{\rho\in L^1(R^3)\cap L^3(R^3):\quad \rho \geq 0,\,\int \rho=1\}$ then
$$
\inf\limits_{\rho\in\mathcal{A}} h(\rho)>-\infty.
$$
Indeed for any $\rho\in \mathcal{A}$ and fixed $R>0$ we have first
\begin{align*}
h(\rho)& =|\rho|_{L^3}-C\int\limits_{|x|\geq R}\frac{\rho}{|x|}-C\int\limits_{|x|\leq R}\frac{\rho}{|x|}\\
& \geq |\rho|_{L^3} -\frac{C}{R}\int\limits_{|x|\geq R}\rho-C\int\limits_{|x|\leq R}\frac{\rho}{|x|}\\
&\geq |\rho|_{L^3} -\frac{C}{R}-C\int\limits_{|x|\leq R}\frac{\rho}{|x|}
\end{align*}
because $\int\limits_{|x|\geq R}\rho\leq |\rho|_{L^{1}(R^3)}=1$. Using now
Young's inequality in the last term with $\rho\in L^{3}(B_R)$ and $|1/|x||_{L^{3/2}(B_R)}=CR$ we get
\begin{align*}
h(\rho)& \geq |\rho|_{L^3} -\frac{C}{R}-C\int\limits_{|x|\leq R}\rho\frac{1}{|x|}\\
& \geq |\rho|_{L^3} -\frac{C}{R}-C|\rho|_{L^{3}(B_R)}|1/|x||_{L^{3/2}(B_R)}\\
& \geq |\rho|_{L^3} -\frac{C}{R} -CR|\rho|_{L^3}.
\end{align*}
Choosing $R>0$ small enough we have therefore
$$
h(\rho)\geq C_1|\rho|_{L^3}-C_2\hspace{3cm}(1)
$$
for some $C_1,C_2>0$ and all $\rho$, and my claim obviously follows.

**Claim 2:** the functional is lower semi-continuous for the weak $L^1\cap L^3$ convergence. Indeed, observe first that $h$ is strongly continuous with respect to strong $L^1\cap L^3$ convergence: the $|\rho|_{L^3}$ term is obviously continuous, and for the second term fix any $R>0$ and argue separately in $B_R$ (use $\frac{1}{|x|}\in L^{3/2}=L^{3'}$ in $B_R$ and strong $L^3$ convergence) and in $R^3\setminus B_R$ (use now $\frac{1}{|x|}\in L^{\infty}$ and strong $L^1$ convergence). Since $h$ is convex my claim is a classical consequence of the Hahn-Banach theorem, see [Brezis, "functional analysis and applications, corollariy 3.9 p. 61]

**Sketch of proof:** By my claim 1 there exists a minimizing sequence $\rho_n\in\mathcal{A}$. By (1) we may assume that $\rho_n\rightharpoonup \rho$ in $L^3(R^3)$, with of course $\rho\geq 0$. If you could prove that $\rho_n\rightharpoonup \rho$ in $L^1(R^3)$ as well we would be done by my claim 2 (the limit $\rho$ would have unit mass since $\int \rho=<\rho,1>_{L^1,L^{\infty}}=\lim\limits_{n\to\infty}<\rho_n,1>_{L^1,L^{\infty}}=\int \rho_n=1$, so that $\rho\in\mathcal{A}$ would be admissible).

Due to the usual lack of $L^1$ reflexivity, obtaining the weak $L^1$ convergence is far from being trivial. By the Dunford-Pettis theorem we know that $\rho_n\rightharpoonup \rho$ iff: (i) $|\rho_n|_{L^1}\leq C$ (OK because they are probability measures), (ii) no concentration (still OK because we have here uniform $L^3$ bounds), and (iii) the familiy $\{\rho_n\}_n\subset L^1$ is uniformly integrable. Uniform integrability means basically that there is no loss of mass at infinity. My feeling is that this should be true, but I can't prove it rigorously. A sketchy argument why this should hold is the following: spreading mass at infinity does not really affect the $L^3$ norm, but forces $\rho$ to "see" $\frac{1}{|x|}$ for larger values of $|x|$. As a consequence spreading mass decreases $\int\frac{\rho}{|x|}$ and increases $h(\rho)$. So the minimizing sequence cannot spread too much mass... Now you need a quantitative result for this claim, maybe trying to control $\int \log(1+ |x|)\rho$ (this is somehow natural because of the $1/r\sim\frac{d}{dr}\log(1+r)$ for large $r$)

One last remark: I imagin you want to derive some PDE for the minimizer $\rho_m$. Because of your non-linear constraint (restriction to probability measures), the differentiation with respect to Euclidean perturbations $\frac{d}{dt}h(\rho+th)$ is not very well adapted. You should try to use instead the optimal transport approach: fix any smooth vector-field $\zeta\in \mathcal{C}^{\infty}_c(R^3,R^3)$ and let $\Phi_{\varepsilon}$ be the associated $\varepsilon$ flow (i-e $d\Phi_{\varepsilon}/d\varepsilon=\zeta(\Phi_\varepsilon)$ and $\Phi_0=\operatorname{Id}$). Looking at push-forward perturbations
$$
\rho_{\varepsilon}:=\{\Phi_\varepsilon\}_{\#}\rho
$$
and computing $\left.\frac{d}{d\varepsilon}h(\rho_\varepsilon)\right|_{\varepsilon=0}=0$ should give you the Euler-Lagrange equation
(see Villani's book "topics in optimal transportation" for a nice introduction). If you had $|\rho|^3_{L^3}$ instead of $|\rho|_{L^3}$ you would end-up with $2\Delta(\rho^3)-c\operatorname{div}(\rho\nabla(1/|x|))=0$, here I don't know what you get. But you should definitly have the $\operatorname{div}(\rho\nabla(1/|x|))$ term.

5more comments