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 the general mass function considering only initial configurations that give a conservative result. We then study the formation and evolution of PBH binaries before recombination by performing N-body simulations. We find that the analytical merger rate estimate based on the tidally perturbed 2-body problem fails when PBHs comprise all dark matter, as most initial binaries are disrupted by the surrounding PBHs. This is 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\).