  We present a two-step method specifically tailored for band structure
calculation of the small-angle moir\'{e}-pattern materials which contain tens
of thousands of atoms in a unit cell. In the first step, the self-consistent
field calculation for ground state is performed with $O(N)$ Krylov subspace
method implemented in OpenMX. Secondly, the crystal momentum dependent Bloch
Hamiltonian and overlap matrix are constructed from the results obtained in the
first step and only a small number of eigenvalues near the Fermi energy are
solved with shift-invert and Lanczos techniques. By systematically tuning two
key parameters, the cutoff radius for electron hopping interaction and the
dimension of Krylov subspace, we obtained the band structures for both rigid
and corrugated twisted bilayer graphene structures at the first magic angle
($\theta=1.08^\circ$) and other three larger ones with satisfied accuracy on
affordable costs. The band structures are in good agreement with those from
tight binding models, continuum models, plane-wave pseudo-potential based
$ab~initio$ calculations, and the experimental observations. This efficient
two-step method is to play a crucial role in other twisted two-dimensional
materials, where the band structures are much more complex than graphene and
the effective model is hard to be constructed.
