In this paper, based on the imaginary time gradient flow model in the density functional theory, a scalar auxiliary variable (SAV) method is developed for the ground state calculation of a given electronic structure system. To handle the orthonormality constraint on those wave functions, two kinds of penalty terms are introduced in designing the modified energy functional in SAV, i.e., one for the norm preserving of each wave function, the other for the orthogonality between each pair of different wave functions. A numerical method consisting of a designed scheme and a linear finite element method is used for the discretization. Theoretically, the desired unconditional decay of modified energy can be obtained from our method, while computationally, both the original energy and modified energy decay behaviors can be observed successfully from a number of numerical experiments. More importantly, numerical results show that the orthonormality among those wave functions can be automatically preserved, without explicitly preserving orthogonalization operations. This implies the potential of our method in large-scale simulations in density functional theory.