Chemical accuracy serves as an important metric for assessing the effectiveness of the numerical method in Kohn–Sham density functional theory. It is found that to achieve chemical accuracy, not only the Kohn–Sham wavefunctions but also the Hartree potential, should be approximated accurately. Under the adaptive finite element framework, this can be implemented by constructing the a posteriori error indicator based on approximations of the aforementioned two quantities. However, this way results in a large amount of computational cost. To reduce the computational cost, we propose a novel solver for the Kohn–Sham equation by using the multi-mesh adaptive technique, in which the Kohn–Sham equation and the Poisson equation are solved in two different meshes on the same computational domain, respectively. With the proposed method, chemical accuracy can be achieved with less computational consumption compared with the adaptive method on a single mesh, as demonstrated in a number of numerical experiments.