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, UO22+, 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.