A method for solving the two-center Dirac equation using a Gaussian-product basis set tailored to the axial symmetry of diatomic systems is presented. For the ground state of H2+ at an internuclear distance R = 2.0 a0, we obtain a relativistic energy with relative precision at the 10−18 level. The properties of the basis set allow for an analytic Fourier transform of the wave function. As a first application we evaluate, for the first time, the zero-potential contribution to the one-loop self-energy for H2+ in Feynman and Coulomb gauges. Although this contribution is not separately gauge invariant, its calculation provides a proof of principle for the use of the present framework in prospective QED calculations.