A common way of representing the long-time dynamics of materials is in terms of a Markov chain that specifies the transition rates for transitions between metastable states. Such chains can either be analyzed directly, used to generate trajectories using kinetic Monte Carlo, or upscaled into mesoscale models such as cluster dynamics. While a number of approaches have been proposed to infer such a representation from direct molecular dynamics (MD) simulations, challenges remain. For example, as chains inferred from a finite amount of MD will in general be incomplete, quantifying their completeness and propagating these uncertainties to observables of interest is extremely desirable. In addition, making the construction of the chain as computationally affordable as possible is paramount. In this work, we simultaneously address these two questions. We first quantify the local completeness of the chain in terms of Bayesian estimators of the yet-unobserved rate, and its global completeness in terms of the residence time of trajectories within the explored subspace. We then systematically reduce the cost of creating the chain by leveraging an accelerated MD method, namely Temperature Accelerated Dynamics. We maximize the increase in residence time against the distribution of states in which additional MD is needed and the temperature at which these are respectively carried out. Using examples of defects that are relevant to the evolution of irradiated materials, we demonstrate that our approach is an efficient, fully automated, and massively-parallel scheme to efficiently explore the long-time behavior of materials.