The diffusion in a binary mixture of species coadsorbed on a multisite square lattice is investigated by a theoretical approach (Chumak and Tarasenko, 1980) based on the non-equilibrium statistical operator method proposed by Zubarev (1961). The investigated lattice gas system is rather complex. There is a mixture of two types of particles (atoms, and/or molecules) adsorbed on a multisite lattice which consists of the three types of adsorption sites. This lattice can be subdivided onto three homogeneous square sublattices composed of the sites of the same type. As the binding energies for the molecules of distinct types adsorbed on the nonidentical sublattices are different, there are six average occupancies (coverages) describing the molecule distribution over the sublattices. A system of the balance equations, which controls the exchange of the molecules between the sublattices on the atomistic level is reduced to the diffusion equations describing the evolution of small hydrodynamic fluctuations of these coverages on the macroscopic level. The diffusion equations are written in the Onsager representation, when the driving forces are gradients of the chemical potentials and in the Fickian representation, when the driving forces are the gradients of coverages. The derivation of these equations results in the analytical expressions for the Fickian diffusivities and Onsager phenomenological coefficients. Despite the complex process of derivation, the final results are simple. The kinetic coefficients, which describe the non-equilibrium state of the molecule system, are expressed via the thermodynamic quantities calculated in the equilibrium state of the system: the chemical potentials (or activities), pair correlation functions and the isothermal susceptibilities. These expressions are derived with account of the lateral interactions between the nearest neighbor molecules. The matrix of the Fickian diffusivities has the ordinary classical form. The diffusivity is a non-symmetric 2 × 2 matrix. The Onsager coefficients are represented by a diagonal 2 × 2 matrix, where the cross terms are absent.