A new method for calculating the ground state electron density of interacting molecules is presented. The supermolecule electron density is obtained using an iterative procedure. At each step the electron density of one molecule is calculated using previously introduced Kohn-Sham equations with constrained electron density. These equations contain terms representing the coupling between constrained and non-constrained electron densities. The coupling terms also involve a new functional, namely the non-additive kinetic energy functional that is not present in the original Kohn-Sham method. Its first-principles analytical form in not yet known. We examine the analytical form of this functional derived from Thomas-Fermi theory. The electron density obtained is compared with that calculated using the original Kohn-Sham method applied to the supermolecule. Good agreement has been found for a broad range of electron density overlaps.