Let $X$ be a metric space with doubling measure, and $L$ be a non-negative, self-adjoint operator satisfying Davies-Gaffney bounds on $L^2(X)$. In this article the authors present a theory of Hardy and BMO spaces associated to $L$, including an atomic (or molecular) decomposition, square function characterization, and duality of Hardy and BMO spaces. Further specializing to the case that $L$ is a Schrodinger operator on $mathbb{R}^n$ with a non-negative, locally integrable potential, the authors establish additional characterizations of such Hardy spaces in terms of maximal functions. Finally, they define Hardy spaces $H^p_L(X)$ for $p>1$, which may or may not coincide with the space $L^p(X)$, and show that they interpolate with $H^1_L(X)$ spaces by the complex method.