  We present for the first time an efficient iterative method to directly solve
the four-component Dirac-Kohn-Sham (DKS) density functional theory. Due to the
existence of the negative energy continuum in the DKS operator, the existing
iterative techniques for solving the Kohn-Sham systems cannot be efficiently
applied to solve the DKS systems. The key component of our method is a novel
filtering step (F) which acts as a preconditioner in the framework of the
locally optimal block preconditioned conjugate gradient (LOBPCG) method. The
resulting method, dubbed the LOBPCG-F method, is able to compute the desired
eigenvalues and eigenvectors in the positive energy band without computing any
state in the negative energy band. The LOBPCG-F method introduces mild extra
cost compared to the standard LOBPCG method and can be easily implemented. We
demonstrate our method in the pseudopotential framework with a planewave basis
set which naturally satisfies the kinetic balance prescription. Numerical
results for Pt$_{2}$, Au$_{2}$, TlF, and Bi$_{2}$Se$_{3}$ indicate that the
LOBPCG-F method is a robust and efficient method for investigating the
relativistic effect in systems containing heavy elements.
