The relativistic Dirac equation for a bound electron in the field of two fixed positive charges is revisited. In contrast to the one center case this three dimensional equation is separable only partially around the azimuthal angle φ, because of the commutation of the Dirac Hamiltonian only with the z component of the total angular momentum Jz. In this work we determine the variational exact solution of this two center problem using a basis constructed by linear combinations of relativistic Slater type spinor wave functions with non integer powers of the radii r1 and r2 on the two centers. We present in some detail the determination of the two center integrations involved. The solutions are obtained by a minimax procedure, that we have developed with a new iterative scheme. We use independent large and small components of the Dirac spinor. This permits us to take control of the spurious solutions, and gives us the possibility to avoid them by the appropriate choice of the wave function parameters. We investigate the behavior of the electron in its 1sσg level of the diatomic homo-nuclear systems A2(2Z-1)+, where A represents the heavy element and Z its atomic number. In the case of heavy ions we study the dependance of the electron energy on the internuclear distance, this gives us an indication for the conditions of atomic collapse, which can induce electron-positron pair production. Our approach has the advantage of needing small basis sets for a relative error of the of 10-7-10-8. It can also be extended easily to the excited level such as the 1sσu level. © 2021 Elsevier B.V.