O(N ) ab initio calculation scheme for large-scale moiré structures

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 the ground state is performed with the $O(N)$ Krylov subspace method implemented in openmx. Second, 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 the Krylov subspace, we obtained the band structures for both rigid and corrugated twisted bilayer graphene structures down to the first magic angle ($\ensuremath{\theta}=1.{08}^{\ensuremath{\circ}}$) with high enough accuracy at affordable costs. The band structures are in good agreement with those from tight-binding models, continuum models, plane-wave pseudopotential based ab initio calculations, and experimental observations. This method is also shown to be efficient in twisted double-bilayer graphene and bilayer ${\mathrm{WSe}}_{2}$. We think this two-step method can play a crucial role in other twisted two-dimensional materials, especially those with much more complex band structure and where the effective model is hard to construct.