  The abundance of primordial black holes (PBHs) in the mass range $0.1 - 10^3
M_\odot$ can potentially be tested by gravitational wave observations due to
the large merger rate of PBH binaries formed in the early universe. To put the
estimates of the latter on a firmer footing, we first derive analytical PBH
merger rate for general PBH mass functions while imposing a minimal initial
comoving distance between the binary and the PBH nearest to it, in order to
pick only initial configurations where the binary would not get disrupted. We
then study the formation and evolution of PBH binaries before recombination by
performing N-body simulations. We find that the analytical estimate based on
the tidally perturbed 2-body system strongly overestimates the present merger
rate when PBHs comprise all dark matter, as most initial binaries are disrupted
by the surrounding PBHs. This is mostly due to the formation of compact N-body
systems at matter-radiation equality. However, if PBHs make up a small fraction
of the dark matter, $f_{\rm PBH} \lesssim 10\%$, these estimates become more
reliable. In that case, the merger rate observed by LIGO imposes the strongest
constraint on the PBH abundance in the mass range $2 - 160 M_\odot$. Finally,
we argue that, even if most initial PBH binaries are perturbed, the present
BH-BH merger rate of binaries formed in the early universe is larger than
$\mathcal{O}(10)\,{\rm Gpc}^{-3} {\rm yr}^{-1}\, f_{\rm PBH}^3$
