A new method is presented to calculate binding energies and eigenfunctions for molecules, using the Dirac–Slater Hamiltonian. A numerical basis set of four component wavefunctions is obtained from atom-like Dirac–Slater wavefunctions. A discrete variational method (DVM) has been applied to generate the binding energies and eigenfunctions for the molecule. Results are given for a series of molecules, including dihydrides H2X (X=O, S, Se, Te), diatomic indium halides InX (X=F, Cl, Br, I), and metal chlorides XCl (X=B, Al, Ga, In, Tl). Comparison is made with results from nonrelativistic calculations using the DVM with numerical Hartree–Fock–Slater-type wavefunctions and with other types of nonrelativistic calculations. In particular, relativistic level shifts and spin–orbit splitting have been analyzed. The theoretical ionization energies are compared with experimental results. Generally a very good agreement is obtained between the experimental and theoretical binding energies for the valence levels, calculated by a transition state procedure. Most of the calculations have been made with a minimal basis set consisting of the occupied atomic orbitals. For H2S calculations are also reported, using an extended basis set. The change in the binding energies for the valence levels is a few tenths of an eV compared with the results using a minimal basis.