A new ATSP2K module is presented for evaluating the electron density function of any multiconfiguration Hartree-Fock or configuration interaction wave function in the non relativistic or relativistic Breit-Pauli approximation. It is first stressed that the density function is not a priori spherically symmetric in the general open shell case. Ways of building it as a spherical symmetric function are discussed, from which the radial electron density function emerges. This function is written in second quantized coupled tensorial form for exploring the atomic spherical symmetry. The calculation of its expectation value is performed using the angular momentum theory in orbital, spin, and quasispin spaces, adopting a generalized graphical technique. The natural orbitals are evaluated from the diagonalization of the density matrix.