An efficient approach to calculate approximate pure-state and transition reduced density matrices in the framework of the multireference relativistic Fock-space coupled cluster (FS CC) theory is proposed. The method is based on the effective operator formalism and consists of the direct substitution of the FS CC Ansatz for a wave operator into the effective operator expression with the subsequent truncation of expansion at the terms quadratic in cluster amplitudes. The final density matrix is defined by active-space density matrices of different ranks ‘dressed’ with contributions from cluster operators. The method gives a connected expression for pure-state density matrices, provided that the intermediate normalisation condition is fulfilled. Moreover, under some additional assumptions, the connectivity can also be ensured for calculated transition property matrix elements and natural transition spinors. The developed technique allows for fast and accurate calculations of one-particle reduced density matrices for a wide range of electronic states. A pilot application of the new technique to construct averaged atomic natural orbital (ANO) basis sets for fully relativistic electronic structure calculations is presented. © 2024 Informa UK Limited, trading as Taylor & Francis Group.