The sixteen-component, no-pair Dirac--Coulomb--Breit equation, derived from
the Bethe--Salpeter equation, is solved in a variational procedure using
Gaussian-type basis functions for the example of positronium, muonium, hydrogen
atom, and muonic hydrogen. The α fine-structure-constant dependence of
the variational energies, through fitting a function of αn and
αnlnα terms, shows excellent agreement with the relevant
energy expressions of the (perturbative) non-relativistic QED framework, and
thereby, establishes a solid reference for the development of a computational
relativistic QED approach