Your operator is self adjoint on $ L^2 (M , e^{-f}) $:
If your operator is $L$.
Then we have for all fonctions u,v
$\int uLv e^{-f}=\int (<du, dv>+huv) e^{-f}$

Hence , the equation $ Lu=0$ has a positive solution
if and only if 0 is the lowest eigenvalue of $L $.

Let's me elaborate on this

This is a consequence of the so called Barta's lemma :
If $u>0$
then a integration by parts formula implies that for every
v:
$$\int_M [|d(uv)|^2+h(uv)^2]e^{-f}= \int_M (Lu)uv^2+|dv|^2 u^2e^{-f}$$
Hence if there is a positive function $u$ such that
$Lu\ge 0$, we get that the spectrum of $L$ is non negative.
(This is the content of Barta's lemma: J. Barta, Sur la vibration fondamentale d'une membrane, {\em C. R. Acad. Sci. Paris} \textbf{204} (1937),
472--473.)

This Barta's lemma implies that if the equation $Lu=0$ has a positive solution, then $0$ is the lowest eigenvalue of $L$.

The other implication follows from the quadratic form :
$$Q(v)=\int_M [|dv|^2+hv^2]e^{-f}= \int_M (Lv)v e^{-f}$$
The lowest eigenvalue , called it $\lambda_0$, of $L$ has a variationnal formulation :
$$\lambda_0 =\inf\{ Q(v), v\in W^{1,2}, \int_M v^2 e^{-f}=1\}$$

And the eigenspace associated to $\lambda_0$ is precisely
$$\{v\in C^\infty(M), Q(v)=\lambda_0\int_M v^2 e^{-f}\}$$
(the fact that these eigenfunctions are smooth is a consequence of Elliptic regularity).

But if $v\in W^{1,2}$ then $|v|\in W^{1,2}$ and $Q(|v|)\le Q(v)$ , $\int_M |v|^2 e^{-f}=\int_M v^2 e^{-f}$, hence if $u$ is an eigenvalue of $L$ associated to $\lambda_0$ then $|u|$ is also an eigenvalue and by linearity
$u+|u|$ and $|u|-u$ also. But the unique continuation properties implies that if $|u|\not =u$ , then the function $|u|+u$ vanishes on the open subset
where $u$ is non positive hence we must have $u+|u|=0$, hence $u$ is non negative. So that there is always a non negative eigenfunction associated to $\lambda_0$. The fact that such a non negative eigenfunction is positive is a consequence of the Harnack inequality (see Gilbarg-Trudinger subsection 8.8).

There is another approach that is more abstract (in the context of Dirichlet form see : The Allegretto-Piepenbrink Theorem
for Strongly Local Dirichlet Forms, by D. Lenz, P. Stollmann, I.Veselic publish in {\em Documenta Mathematica} \textbf{14} (2009) 167--189).

For complete non compact manifold, the existence of a positive solution of $Lu=0$ is equivalent to the fact that the spectrum of $L$ is non negative. This is called the Allegretto-Piepenbrink principle and has been adapted to manifold by Fischer-Colbrie and Schoen
(see W.F. Moss and J. Piepenbrink. Positive Solutions of Elliptic Equations. {\em Pacific J. Math.},
\textbf{75}:219--226, 1978.

D. Fischer-Colbrie and R. Schoen. The Structure of Complete Stable Minimal Surfaces in
3-manifolds of Non-negative Scalar Curvature.{\em Comm. Pure Appl. Math.}, \textbf{XXXIII}:199--211,
1980.

and the lemma 3.10 in the beautiful book : S. Pigola, M. Rigoli, and A. Setti. {\em Vanishing and Finiteness Results in Geometric Analysis.}
Birkhäuser, 2008.