We show that all-electron relativistic four-component Dirac-Kohn-Sham (DKS) computations, using G-spinor basis sets and state-of-the-art density fitting algorithms, can be efficiently parallelized and applied to large molecular systems, including large clusters of heavy atoms. The performance of the parallel implementation of the DKS module of the program BERTHA is illustrated and analyzed by some test calculations on several gold clusters up to Au32, showing that calculations with more than 25 000 basis functions (i.e., DKS matrices on the order of 10 GB) are now feasible. As a first application of this novel implementation, we investigate the interaction of the atom Hg with the Au20 cluster.