I think this has to do with what you treat as the dynamic variables. In your first formulation, $\rho$ was given and fixed, and you wanted to calculate $V$ by varying only $V$, not $\rho$. In the second formulation, $V$ is just a quantity derived from $\rho$, and this is the energy suitable for obtaining $\rho$ by varying $\rho$, with $V$ just a shorthand for a certain linear transformation of $\rho$. A third possibility is to consider $V$ as given and find the motion of a charge distribution in this potential -- here, again, there would not be a factor of $\frac{1}{2}$. This is the case, for instance, if we calculate the motion of the electron in the potential of a hydrogen nucleus, or of the Earth in the gravitational potential of the Sun.
Mathematically speaking, we need a factor of $\frac{1}{2}$ in front of quadratic terms (where $V\rho$ is quadratic in $\rho$ if $V$ is considered as a derived quantity derived from the dynamical variable $\rho$) but not in front of linear terms. Physically speaking, the factor of $\frac{1}{2}$ avoids double-counting when the energy is regarded as the interaction energy of a charge distribution with itself.
Now you may ask: But surely there is some well-defined energy and we can't just pick and choose the factors according to expedience? That's true, but consider again the hydrogen atom: If we consider the potential of the nucleus as given, then we're not counting the potential energy of the nucleus in the field of the electron in the integral over $\rho V$. But the entire energy is there and has to be accounted for in the integral, hence the factor $1$. On the other hand, in treating a helium atom, where both electrons are treated as dynamic and their interaction energy enters into the Hamiltonian, we do include a factor of $\frac{1}{2}$ to avoid double-counting the energy in the double integral over $\rho_1\rho_2/r_{12}$.
It is not appropriate to work in the space $H^1_0(U)$ since nonzero boundary conditions are being considered. You will have to assume some regularity of the boundary of $U$.
One version of Green's theorem (see e.g. the appendix in Evans) is that $$- \int_U (\Delta u) v \, dx = \int_U Du \cdot Dv \, dx -\int_{\partial U} \frac{\partial u}{\partial \nu} v \, dS.$$
A weak solution to the problem at hand can be proposed by setting $-\Delta u = f$ in $U$ and $\dfrac{\partial u}{\partial \nu} = -u$ on $\partial U$ so that $$\int_U fv = \int_U Du \cdot Dv + \int_{\partial U} uv \, dS \quad \forall v \in H^1(U)$$
or a bit more precisely
$$\int_U fv = \int_U Du \cdot Dv + \int_{\partial U} \newcommand{\tr}{\mathrm{Tr}\ \! }( \tr u )(\tr v) \, dS \quad \forall v \in H^1(U)$$
where $\tr : H^1(U) \to L^2(\partial U)$ is the trace operator.
An appropriate bilinear form is thus given by $$B[u,v] = \int_U Du \cdot Dv + \int_{\partial U} \newcommand{\tr}{\mathrm{Tr}\ \! }( \tr u )(\tr v) \, dS, \quad u,v \in H^1(U).$$
$B$ is clearly bounded. As far as coercivity goes, it may be helpful to use the Rellich-Kondrachov theorem. I can follow up with a hint if you like.
It remains to show that there is a constant $\alpha > 0$ with the property that $\|u\|_{H^1}^2 \le \alpha B[u,u]$ for all $u \in H^1(U)$. This can be proven by contradition. Otherwise, for every $n \ge \mathbb N$ there would exist $u_n \in H^1(U)$ with the property that $\|u_n\|^2_{H^1} > n B[u_n,u_n]$. For each $n$ define $v_n = \dfrac{u_n}{\|u_n\|_{H^1}}$. Then $v_n \in H^1(U)$, $\|v_n\|_{H^1} = 1$, and $B[v_n,v_n] < \dfrac 1n$ and all $n$.
Here we can invoke Rellich-Kondrachov. Since the family $\{v_n\}$ is bounded in the $H^1$ norm, there is a subsequence $\{v_{n_k}\}$ that converges to a limit $v \in L^2(U)$. However, since $\|Dv_{n_k}\|_{L^2}^2 < \dfrac{1}{n_k}$ it is also true that $Dv_{n_k} \to 0$ in $L^2$. Thus for any $\phi \in C_0^\infty(U)$ you have $$\int_U v D \phi \, dx = \lim_{k \to \infty} \int_U v_{n_k} D \phi \, dx = - \lim_{k \to \infty} \int_U D v_{n_k} \phi \, dx = 0.$$ This means $v \in H^1(U)$ and $D v = 0$, from which you can conclude $v_{n_k} \to v$ in $H^1(U)$. Since $\|v_{n_k}\|_{H^1} = 1$ for all $k$ it follows that $\|v\|_{H^1} = 1$ as well.
Next, since $\|\tr v_{n_k}\|_{L^2(\partial U)}^2 < \dfrac{1}{n_k}$ and the trace operator is bounded there is a constant $C$ for which
$$ \|\tr v\|_{L^2(\partial \Omega)} \le \|\tr v - \tr v_{n_k}\|_{L^2(\partial \Omega)} + \|\tr v_{n_k}\|_{L^2(\partial \Omega)} < \frac{1}{n_k} + C \|v - v_{n_k}\|_{H^1(U)}.$$
Let $k \to \infty$ to find that $\tr v = 0$ in $L^2(\partial U)$.
Can you prove that if $v \in H^1(U)$, $Dv = 0$, and $\tr v = 0$, then $v = 0$? Once you have established that fact you arrive at a contradiction, since $v$ also satisfies $\|v\|_{H^1} = 1$. It follows that $B$ is in fact coercive.
Best Answer
Short answer. If $\Omega$ is unbounded, you still have uniqueness but you might lose existence in $H^1_0(\Omega)$.
Long answer. You need boundedness of $\Omega$ to use the typical argument based on Riesz's representation of linear functionals (that's the argument that Jose mentions in his comment). This is because Poincaré's inequality $c\int_{\Omega}u^2\, dx \le \int_\Omega \lvert \nabla u\rvert^2\, dx$ is not available on unbounded domains.
However, you still have uniqueness in $H^1_0(\Omega)$, regardless of whether $\Omega$ is bounded or not. If $u_1$ and $u_2$ are two $H^1_0(\Omega)$ weak solutions to $$\tag{$\star$} \begin{cases}-\Delta u = f, & \Omega \\ u=0, & \partial \Omega\end{cases} $$ then the difference $v=u_1-u_2$ satisfies $$ \int_\Omega \lvert \nabla v\rvert^2\, dx =0$$ and so it is a constant. (I am here assuming that domain means connected open set). The only constant function in $H^1_0(\Omega)$ is the null function. Therefore $v=0$.
Remark 1. If you do not impose the condition $u\in H^1_0(\Omega)$, but rather look for classical solutions to $(\star)$, then you do lose uniqueness in unbounded domains. For example, in $\Omega=\mathbb{R}^n$ you can add any entire harmonic function to a solution and obtain another solution. The point is that $(\star)$ is underdetermined in unbounded domains, because no boundary condition is imposed at infinity. By insisting that the solution be in $H^1_0(\Omega)$ one forces some decay at infinity and thus recovers uniqueness.
Remark 2. The existence of solutions to $(\star)$ in unbounded domains is another story. Under the sole assumption that $f\in L^2(\Omega)$ it might happen that no $H^1_0(\Omega)$ weak solution exist. For a concrete example, fix $n= 3$ and $\Omega=\mathbb{R}^3$. Let $$ u(x)=\left(\frac{1}{1+\lvert x \rvert^2}\right)^\frac{1}{2}$$ and $$ f(x)=3\left(\frac{1}{1+\lvert x \rvert^2}\right)^\frac{5}{2}.$$ (The function $u$ appears in the elliptic theory of PDEs. It is called standard bubble). You can check that $$-\Delta u = f$$ but $u\notin L^2(\mathbb{R}^3)$ (hence $u\notin H^1(\mathbb{R}^3)$) even if $f\in L^2(\mathbb{R}^3)$. And $u$ is the unique solution to $(\star)$ with the given source term $f$, at least in the class of tempered distributions. In particular, no $H^1(\mathbb{R}^3)$ solutions to $(\star)$ with the given source term $f$ exist.
To see things clearly in $\Omega=\mathbb{R}^n$ you can solve $(\star)$ by means of the Fourier transform (or by convolution against the fundamental solution of the Laplace operator). Then the solution is formally written at Fourier side as $$\hat{u}=\frac{\hat{f}}{\lvert\xi\rvert^2}, $$ and this is the unique solution to $(\star)$ in the class of tempered distributions, as announced previously. Now it happens that $u\in H^1_0(\mathbb{R}^n)=H^1(\mathbb{R}^n)$ if and only if $$ \frac{(1+\lvert\xi\rvert^2)^\frac{1}{2}}{\lvert\xi\rvert^2} \hat{f}\in L^2(\mathbb{R}^n).$$ This is not guaranteed by the condition $f\in L^2(\mathbb{R}^n)$ alone: to have this, $\hat{f}$ must vanish at $\xi=0$ to second order at least.