  We present the first full-potential method that solves the fully relativistic
4-component Dirac-Kohn-Sham equation for materials in the solid state within
the framework of atom-centered Gaussian-type orbitals (GTOs). Our GTO-based
method treats one-, two-, and three-dimensional periodic systems on an equal
footing, and allows for a seamless transition to the methodology commonly used
in studies of molecules with heavy elements. The scalar relativistic effects
and spin-orbit interaction are handled variationally. The full description of
the electron-nuclear potential in the core region of heavy nuclei is
straightforward due to the local nature of the GTOs and does not pose any
computational difficulties. We show how the time-reversal symmetry and a
quaternion algebra-based formalism can be exploited to significantly reduce the
increased methodological complexity and computational cost associated with
multiple wave-function components coupled by the spin-orbit interaction. We
provide a description of how to employ the matrix form of the multipole
expansion and an iterative renormalization procedure to evaluate the
conditionally convergent lattice sums arising for periodic systems. We
investigate the problem of inverse variational collapse that arises if the
Dirac operator containing a repulsive periodic potential is expressed in a
basis containing diffuse functions, and suggest a possible solution. We
demonstrate the validity of the method on silver halide crystals with large
relativistic effects, and two-dimensional honeycomb structures (silicene and
germanene) exhibiting the spin-orbit-driven quantum spin Hall effect. Our
results are well-converged with respect to the basis limit using standard bases
developed for molecular calculations, and indicate that the rule of removing
basis functions with small exponents should not be applied when transferring
the molecular basis to solids.
