Electronic band structure, ferroelectric and magnetic properties of Rhombohedral (hexagonal axis) phase of multiferroic BiFeO3 and RE (RE = La, Nd, Sm and Gd) doped compounds have been investigated using TB-mBJ semi local (Tran-Blaha modified Becke-Johnson) potential approximation method using WIEN2k code. Studies were also undertaken to understand the influence of spin-orbit coupling (SOC) corrections on all the properties of these materials. For this purpose, the total energies were calculated as a function of reduced volumes and the data were fitted to Brich Murnaghan equation. The estimated ground state parameters are found to be comparable with those of experimental ones. The semiconducting band gap values of all the samples were also obtained using both the TB-mBJ and TB-mBJ + SOC methods. One may understand from the density of states that the valence band consists of Fe-d and O-p states and the conduction band is having Fe-d and Bi-p states. From the charge density studies one may see that stereo chemically active lone-pair electrons are present at Bi sites and are responsible for the displacements of Bi atoms from the centro-symmetric to the non-centrosymmetric structure resulting in the exhibition of ferroelectricity. Further, it has been concluded that "lone pair" might have formed due to hybridization of 6s and 6p atomic orbitals with 6s(2) electrons filling one of the resulting orbitals in Bi. The calculated magnetic moments at the iron sites are not integer values probably due to hybridization of Fe electrons with the neighboring O ions.