An implementation of a four-component density corrected approach for calculations of nuclear electric field gradients (EFGs) in molecules based on the two-component relativistic zeroth-order regular approximation (ZORA) is reported. The program module, which is part of the NWChem package, allows for scalar and spin-orbit relativistic computations of EFGs. Benchmark density functional calculations are reported for a large set of main group diatomic molecules, a set of Cu and Au diatomics, several Ru and Nb complexes, the free uranyl ion, and two uranyl carbonate complexes. Data obtained from nonhybrid as well as fixed and range-separated hybrid functionals are compared. To allow for a chemically intuitive interpretation of the results, a breakdown of the EFGs of selected systems in terms of localized molecular orbitals is given. For CuF, CuCl, AuCl, UO2(2+), and a uranyl carbonate complex, the localized orbital decomposition demonstrates in particular the role of the valence metal d and f shells, respectively, and leads to rather compact analyses. For f orbitals, a Townes-Dailey-like model is set up to assist the analysis.