Let's consider this boundary-value problem:
$$\begin{cases} -\Delta V = \rho & \rm{in}\ \Omega \\ V=0 & \rm{on}\ \partial \Omega \end{cases}.$$
We know that this problem has a variational formulation: its solutions are critical points of the energy functional
$$I[u]=\int_{\Omega} \left( \frac{1}{2}\lvert \nabla u \rvert^2 - u\rho \right)\, dx,$$
(cfr. Evans Partial differential equations 2nd ed., §2.2.5).
Now let us look at one (of the many) possible physical interpretations of this problem, the electrostatic one. If units are chosen so that $\varepsilon_0=1$, this equation describes the electric potential $V$ in a region $\Omega$ filled with a charge distribution $\rho$ and whose boundary is grounded.
It would be nice if the energy functional introduced above coincided with the total energy of this physical system. Unfortunately, it seems to me not to be so. In fact, as I read in Feynman's Lectures on physics, vol.II, §8-5, the total energy of this system consists of two additive parts: one is
$$\frac{1}{2}\int_{\Omega} \lvert \nabla V \rvert^2\, dx$$
and is due to the electric field; the other is
$$\frac{1}{2}\int_{\Omega} V\rho\, dx$$
and is due to the charge distribution.
This induces me to think that the energy functional should be
$$U[u]=\int_{\Omega} \left(\frac{1}{2}\lvert \nabla u \rvert^2 + \frac{1}{2}u\rho \right)\, dx,$$ instead of the correct $$I[u]=\int_{\Omega} \left( \frac{1}{2}\lvert \nabla u \rvert^2 - u\rho \right)\, dx.$$ Why am I wrong?
In the Lagrangian formulation of classical mechanics, the functional that we extremise to derive the Euler-Lagrange equation of motion is the Lagrangian, defined as
$$L[u] = T[u] - V[u]$$
where $T[u]$ is the kinetic part and $V[u]$ is the potential part. Contrast this with the Hamiltonian $H=T+V$ from which we derive Hamilton's equations of motion (these are equivalent to the Euler-Lagrange equations, and the two are related by Legendre transform). That should explain the minus sign discrepancy.
I am not sure about the factor of 1/2. Is it possible that we are simply defining $\rho_{\mathrm{Feynman}} = 2\rho_{\mathrm{Evans}}$?