Basis of $H_0^1$ of solutions to the time-independent Schrödinger equation Let $\Omega\subset \mathbb{R}^N$ be a bounded $C^1$ domain, $q\in L^{\infty}(\Omega;\mathbb{R}_+)$. A function $u\in H_0^1(\Omega)$ is called a weak solution to 
$$(*)\begin{cases}-\Delta u + qu = \lambda u &\text{in }\Omega\\ 
u=0& \text{in }\partial \Omega
\end{cases} $$
for some $\lambda\in \mathbb{R}$, iff it satisfies 
$$\int \nabla u \nabla h= \int (\lambda-q)uh,\qquad \forall h\in H_0^1(\Omega) $$
I am looking for a Hilbert basis $\left\{v_n\right\}$ of $H_0^1$ such that each $v_n$ is an eigenvector of the linear operator $T:L^2(\Omega)\to L^2(\Omega)$, where $T$ maps each $f\in L^2(\Omega)$ to the unique weak solution $u\in H_0^1(\Omega)$ of the problem 
$$(**)\begin{cases}-\Delta u + qu = f& \text{in }\Omega\\ 
u=0& \text{in }\partial \Omega\end{cases} $$
I know that $T$ is compact and positive so by the Hilbert-Schmidt theorem, I have an orthonormal basis $\left\{\varphi_n\right\}$ of eigenvectors of $T$ for $L^2(\Omega)$, with positive eigenvalues $\left\{\lambda_n^{-1}\right\}$ with $\lambda_n\to +\infty$.  Moreover, each $\varphi_n$ is in $H_0^1(\Omega)$, and it satisfies 
$$\int \nabla u \nabla h=\int (\lambda_n-q)uh ,\qquad \forall h\in H_0^1$$
i.e. it is a weak solution to (*) with $\lambda=\lambda_n$.
The idea is to prove that $\left\{\varphi_n\right\}$ is also a basis for $H_0^1$, the issue is that the inner product is different, and poses a more restrictive condition to convergence. If $q=0$, then the solution is obvious, because the above equality says that 
$$\left \langle \nabla \varphi_n,\nabla \varphi_m\right\rangle=\lambda_n \left \langle \varphi_n,\varphi_m\right\rangle =\lambda_n \delta_{nm} $$
Therefore $\left\{\varphi_n\right\}$ is also a (non-orthonormal) basis for $H_0^1$. But how does one handle the case of a generic $q\in L^{\infty}(\Omega;\mathbb{R}_+)$, which messes up all the inner products? 
 A: There are two key points that you are overlooking in your analysis of the case $q=0$ that actually solve the problem in general.  
The issue first is that the set $\{\varphi_n\}$ is NOT an orthonormal basis, but just an orthogonal one.  Indeed, we have that 
$$
\int \nabla \varphi_n \cdot \nabla \varphi_m + q \varphi_n \varphi_m = \lambda_n \int \varphi_n \varphi_m = \lambda_n \delta_{nm}.
$$
In particular, when $n=m$ we have that 
$$
\int |\nabla \varphi_n|^2  + q |\varphi_n|^2 = \lambda_n,
$$
which in general means that the functions are not normalized.  To fix this we simply divide by $\sqrt{\lambda_n}$.  
Now we arrive at the second issue: the choice of inner-product.  Plugging in above shows that $\{\varphi_n /\sqrt{\lambda_n}\}$ is indeed an orthonormal basis of $H^1_0$ with respect to the inner-product 
$$
(u,v) = \int \nabla u \cdot \nabla v + q uv.
$$
Since you are assuming that $q \ge 0$, the Poincare inequality tells us that this inner-product gives us a norm that is equivalent to the usual one on $H^1_0$.  The take-away is that if we want orthonormality of the eigenfunctions, then we are forced to view this in the context of the inner-product related to the PDE.  As we change $q$, we change the PDE and the inner-product, but this always generates a norm equivalent to the usual one, so we don't change the underlying topology.
