  We have developed a linear scaling algorithm for calculating
maximally-localized Wannier functions (MLWFs) using atomic orbital basis. An
O(N) ground state calculation is carried out to get the density matrix (DM).
Through a projection of the DM onto atomic orbitals and a subsequent O(N)
orthogonalization, we obtain initial orthogonal localized orbitals. These
orbitals can be maximally localized in linear scaling by simple Jacobi sweeps.
Our O(N) method is validated by applying it to water molecule and wurtzite ZnO.
The linear scaling behavior of the new method is demonstrated by computing the
MLWFs of boron nitride nanotubes.
