We derive a density matrix (DM) theory for quantum cascade lasers (QCLs) that describes the influence of scattering on coherences through a generalized scattering superoperator. The theory enables quantitative modeling of QCLs, including localization and tunneling effects, using the well-defined energy eigenstates rather than the ad hoc localized basis states required by most previous DM models. Our microscopic approach to scattering also eliminates the need for phenomenological transition or dephasing rates. We discuss the physical interpretation and numerical implementation of the theory, presenting sets of both energy-resolved and thermally averaged equations which can be used for detailed or compact device modeling. We illustrate the theorys applications by simulating a high performance resonant-phonon terahertz (THz) QCL design which cannot be easily or accurately modeled using conventional DM methods. We show that the theorys inclusion of coherences is crucial for describing localization and tunneling effects consistent with experiment.