In the framework of Thomas-Fermi approximation, we study systematically the EOSs and microscopic structures of neutron star matter in a vast density range with $n_\mathrm{b}\approx 10^{-10}$-2 $\mathrm{fm}^{-3}$, where various covariant density functionals are adopted, i.e., those with nonlinear self couplings (NL3, PK1, TM1, GM1, MTVTC) and density-dependent couplings (DD-LZ1, DDME-X, PKDD, DD-ME2, DD2, TW99). It is found that the EOSs generally coincide with each other at $n_\mathrm{b}\lesssim 10^{-4}$ fm${}^{-3}$ and 0.1 fm${}^{-3}\lesssim n_\mathrm{b} \lesssim 0.3$ fm${}^{-3}$, while in other density regions they are sensitive to the effective interactions between nucleons. By adopting functionals with larger slope of symmetry energy $L$, the curvature parameter $K_\mathrm{sym}$ and neutron drip density generally increase, while the droplet size, proton number of nucleus, core-crust transition density, and onset density of non-spherical nuclei decrease. All functionals predict neutron stars with maximum masses exceeding the two-solar-mass limit, while those of DD2, DD-LZ1, DD-ME2, and DDME-X predict optimum neutron star radii according to the observational constraints. Nevertheless, the corresponding skewness coefficients $J$ are much lager than expected, while only the functionals MTVTC and TW99 meet the start-of-art constraints on $J$. More accurate measurements on the radius of PSR J0740+6620 and the maximum mass of neutron stars are thus essential to identify the functional that satisfies all constraints from nuclear physics and astrophysical observations. Approximate linear correlations between neutron stars' radii at $M=1.4 M_{\odot}$ and $2 M_{\odot}$, the slope $L$ and curvature parameter $K_\mathrm{sym}$ of symmetry energy are observed as well, which is mainly attributed to the curvature-slope correlations in the functionals adopted here. The results presented here are applicable for the investigations on the structures and evolutions of compact stars in a unified manner.