We present a simple and robust implementation of the phase field fracture method in Abaqus. Unlike previous works, only a user material (UMAT) subroutine is used. This is achieved by exploiting the analogy between the phase field balance equation and heat transfer, which avoids the need for a user element mesh and enables taking advantage of Abaqus' in-built features. A unified theoretical framework and its implementation are presented, suitable for any arbitrary choice of crack density function and fracture driving force. Specifically, the framework is exemplified with the so-called AT1, AT2 and phase field-cohesive zone models (PF-CZM). Both staggered and monolithic solution schemes are handled. We demonstrate the potential and robustness of this new implementation by addressing several paradigmatic 2D and 3D boundary value problems. The numerical examples show how the current implementation can be used to reproduce numerical and experimental results from the literature, and efficiently capture advanced features such as complex crack trajectories, crack nucleation from arbitrary sites and contact problems. The code developed can be downloaded from www.empaneda.com/codes.