Linear scaling density functional theory with Gaussian orbitals and periodic boundary conditions

We report methodological and computational details of our Kohn-Sham density functional method with Gaussian orbitals for systems with periodic boundary conditions (PBC). When solving iterative self-consistent field (SCF) equations of density functional theory (DFT), the most computationally demanding tasks are Kohn-Sham (or Fock) matrix formation and the density matrix update step. The former requires evaluation of the Coulomb interactions and the exchange-correlation quadrature, and in our code both of them are computed via $O$(N) techniques. An $O$(N) approach for the Coulomb problem in electronic structure calculations with PBC is developed here and is based on the direct space fast multipole method (FMM). The FMM achieves not only linear scaling of computational time with system size but also high accuracy, which is pivotal for avoiding numerical instabilities that have previously plagued calculations with large bases, especially those containing diffuse functions. The density matrix update step is carried out via the conventional $O$(N3) diagonalization of the Fock matrix, which for systems with less than ≈3000 basis functions is cheaper than the recently developed $O$(N) algorithms. In addition to evaluating energy, our code also computes analytic energy gradients with respect to atomic positions and cell dimensions (forces). Combining the latter with the developed in this work redundant internal coordinate algorithm for optimization of periodic systems, it becomes possible to optimize geometries of periodic structures with great efficiency and accuracy. We demonstrate the capabilities of our method with benchmark calculations on polyacetylene, poly(p-phenylenevinylene) (PPV), and a series of carbon and boron-nitride single wall nanotubes employing basis sets of double zeta plus polarization quality, in conjunction with generalized gradient approximation and kinetic energy density dependent functionals. We also present vibrational frequencies for PPV obtained from finite differences of forces. The largest calculation reported in this work contains 244 atoms and 1344 contracted Gaussians in the unit cell.