During collisions of heavy ions with heavy targets below the Coulomb barrier, adiabatic molecular orbitals are formed for the inner electrons. Deviations from adiabaticity lead to coupling between various states and can be treated by time-dependent perturbation theory. For high charges ( Z1+Z2 ≧ 60) the molecular electrons are highly relativistic. Therefore, the Dirac equation has to be used to obtain the energies and wave functions. The Dirac Hamiltonian is transformed into the intrinsic rotating coordinate system where prolate spheroidal coordinates are introduced. A set of basis functions is proposed which allows the evaluation of all matrix elements of the Dirac Hamiltonian analytically. The resulting matrix is diagonalized numerically. The finite nuclear charge distribution is also taken into account. Results are presented and discussed for various characteristic systems, e. g. Br-Br, Ni-Ni, I-I, Br-Zr, I-Au, U -U, etc.