Those equations are fundamental models in the simulation of plasma physics. The total energy is an important physical quantity that is conserved by those models. Our methods are the first Eulerian solver that can preserve fully discrete total energy conservation. The main features of our methods include energy-conservative temporal and spatial discretization. In particular, an energy-conserving operator splitting is proposed to enable efficient calculation of fully implicit methods. We validate our schemes by rigorous derivations and benchmark numerical examples.