We present a numerical method that evolves a two-temperature, magnetized, radiative, accretion flow around a black hole, within the framework of general relativistic radiation magnetohydrodynamics. As implemented in the code KORAL, the gas consists of two sub-components ions and electrons - which share the same dynamics but experience independent, relativistically consistent, thermodynamical evolution. The electrons and ions are heated independently according to a prescription from the literature for magnetohydrodynamical turbulent dissipation. Energy exchange between the particle species via Coulomb collisions is included. In addition, electrons gain and lose energy and momentum by absorbing and emitting synchrotron and bremsstrahlung radiation and through Compton scattering. All evolution equations are handled within a fully covariant framework in the relativistic fixed-metric space-time of the black hole. Numerical results are presented for five models of low-luminosity black hole accretion. In the case of a model with a mass accretion rate M similar to 4 x 10(-8) M-Edd, we find that radiation has a negligible effect on either the dynamics or the thermodynamics of the accreting gas. In contrast, a model with a larger M similar to 4 x 10(-4) M-Edd behaves very differently. The accreting gas is much cooler and the flow is geometrically less thick, though it is not quite a thin accretion disc.