Abstract. A four-component density functional program package (Beijing Density Functional), suitable for the calculation of total-energy-related chemical properties of systems containing heavy atoms, was developed. The code is based on modern sophisticated exchange-correlation functionals and was applied to calculate the spectroscopic constants of the lanthanide diatomic molecules of EuO, EuS, YbO and YbS. It is suggested that the experimental bond lengths for EuS and YbS, derived from empirical interpolations, need to be revised. Relativistic effects on the electronic structure are discussed and compared with results from previous work. The involvement of the lanthanide valence orbitals in chemical bonding is investigated with a newly developed population and bonding analysis approach.