I stumbled upon the following inequality where $f \in C_c^{\infty}(\mathbb{R}^3)$:
$$ \left( \int \int_{\mathbb{R}^6} \frac{\overline{f(x)} f(y)}{\vert x-y \vert}dxdy \right) \Vert\nabla f\Vert_2^2 \geq 4 \pi \Vert f\Vert_2^4 $$
where $\overline{f(x)}$ denotes complex conjugation. Unfortunately, I have no idea how to prove it. I don't even know why the left-hand side should be non-negative.
The symmetry of the left-hand side and the $4 \pi$ on the right suggests some kind of transformation into a polar coordinates or use of the co-area formula.
I'm also curious about generalizations to $\mathbb{R}^d$, i.e. if something like
$$ \left( \int \int_{\mathbb{R}^{2d}} \frac{\overline{f(x)} f(y)}{\vert x-y \vert}dxdy \right) \Vert\nabla f\Vert_2^2 \geq \omega_n \Vert f\Vert_2^4 $$
holds for $f \in C_c^{\infty}(\mathbb{R}^d)$, where $\omega_n$ is the surface area of the $d$-dimensional sphere.
It is well known that the Fourier transform of radial functions on $\mathbb{R}^n$ satisfy, for $\alpha \in (0,n)$,
$$ \widehat{\frac{1}{|x|^{n - \alpha}}} = \frac{c}{|\xi|^\alpha} $$
where both sides are understood in the sense of tempered distributions and $c$ is some constant depending on $\alpha$ and $n$ (see, e.g., Theorem 4.1 in Chapter IV of Stein and Weiss, Introduction to Fourier Analysis on Euclidean Spaces). This means that
$$ \int_{\mathbb{R}^3} \frac{f(y)}{|x-y|} ~\mathrm{d}y = [\frac{1}{|\cdot|} * f] (x) $$
has Fourier transform
$$ c \hat{f}(\xi) / |\xi|^2 $$
So by Plancherel's identity you have
$$ \iint_{\mathbb{R}^6} \frac{\bar{f}(x) f(y)}{|x-y|} ~\mathrm{d}y ~\mathrm{d}x = c \int_{\mathbb{R}^3} \frac{|\hat{f}(\xi)|^2}{|\xi|^2} ~\mathrm{d}\xi $$
is a positive real number.
The inequality you are looking at is really a consequence of
$$ \int |\hat{f}(\xi)|^2 ~\mathrm{d}\xi = \int |\xi|^2 \cdot |\xi|^{-2} |\hat{f}(\xi)|^2 ~\mathrm{d}\xi \leq \left( \int |\xi|^2 |\hat{f}(\xi)|^2 ~\mathrm{d}\xi \right)^{\frac12} \cdot \left( \int |\xi|^{-2} |\hat{f}(\xi)|^2 ~\mathrm{d}\xi \right)^{\frac12} $$
by Cauchy-Schwarz.
The correct higher dimensional generalizations are
$$ \left( \iint_{\mathbb{R}^d \times \mathbb{R}^d} \frac{\bar{f}(x) f(y)}{|x-y|^{\color{red}{d - 2}}} ~\mathrm{d}y ~\mathrm{d}x \right) \cdot \| \nabla f\|^2_2 \geq c \| f\|_2^4 $$
where you are probably correct about $c = \omega_n$, but I don't have the time to check the values of the gamma function to make sure.
Side remark: the Fourier transform property described above for the operation of convolving a function against $1 / |x|^\alpha$ for $\alpha \in (0,n)$ earns these operations the name of "fractional integral operators". They are the "inverse" operations to taking (possibly fractional) powers of the Laplacian. So another way to write your inequality is the less-surprising-looking $$ \| |\nabla|^{-1} f \|_{2} \|\nabla f\|_{2} \geq \|f\|_2^2 $$