  We explore the possibility that Dark Matter (DM) may be explained by a
non-uniform background of approximately stellar-mass clusters of Primordial
Black Holes (PBHs), by simulating the evolution them from recombination to the
present with over 5000 realisations using a Newtonian $ N $-body code. We
compute the cluster rate of evaporation, and extract the binary and merged
sub-populations along with their parent and merger tree histories, lifetimes
and formation rates; the dynamical and orbital parameter profiles, the degree
of mass segregation and dynamical friction, and power spectrum of close
encounters. Overall, we find that PBHs can constitute a viable DM candidate,
and that their clustering presents a rich phenomenology throughout the history
of the Universe. We show that binary systems constitute about 9.5\% of all PBHs
at present, with mass ratios of $ \bar{q}_{\rm B} = 0.154 $, and total masses
of $ \bar{m}_{\rm T,\,B} = 303\,M_\odot$. Merged PBHs are rare, about 0.0023\%
of all PBHs at present, with mass ratios of $ \bar{q}_{\rm B}= 0.965 $ with
total and chirp masses of $ \bar{m}_{\rm T,\,B}= 1670\,M_\odot$ and $
\bar{m}_{c,{\rm M}} = 642\,M_\odot $ respectively. We find that cluster puffing
up and evaporation leads to bubbles of these PBHs of order 1 kpc containing at
present times about 36\% of objects and mass, with hundred pc sized cores. We
also find that these PBH sub-haloes are distributed in wider PBH haloes of
order hundreds of kpc, containing about 63\% of objects and mass, coinciding
with the sizes of galactic halos. We find at last high rates of close
encounters of massive Black Holes ($ M \sim 1000\,M_\odot$), with $
\Gamma^{\mathrm{S}} = (1.2^{+5.9}_{-0.9}) \times 10^{7} \mathrm{yr^{-1}
Gpc^{-3}}$ and mergers with $\Gamma^{\mathrm{M}} = 1337 \pm 41 \mathrm{yr^{-1}
Gpc^{-3}} $.
