Abstract. 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 we develop a theory of Hardy and BMO spaces associated to L, including an atomic (or molecular) decomposition, square function characterization, duality of Hardy and BMO spaces. Further specializing to the case that L is a Schrödinger operator on R n with a non-negative, locally integrable potential, we establish addition characterizations of such Hardy space in terms of maximal functions. Finally, we 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.