High order computational approach for the electronic structure calculation is a long-standing research area. Among all the discretization methods, spectral element method stands out for its flexibility to capture the singularity and high accuracy from the spectral perspective. In this work, we propose a tetrahedral spectral element method for full-potential electronic structure calculation. In the method, a C space is constructed for the numerical solution, based on a modification of basis functions proposed in [L. Jia et al. (2022) [48]] with the aid of generalized Koornwider polynomials, and a special designed tetrahedral mesh. Moreover, the interpolation and visualization are developed to provide numerical insights of our method. To demonstrate its potential towards practical simulations, we apply our method for calculating ground states of given electronic structures using the Kohn-Sham model. The method consists of a combination of an imaginary time propagation method and a self-consistent field iteration. For accelerating the simulation, a p-multigrid solver is developed which takes the advantage of the high order property of our method. A number of numerical experiments validate the effectiveness of our method, i.e., the expected spectral accuracy can be observed from results, and convincing ground states of various atoms and molecules are obtained.