The LIGO/Virgo collaboration have by now detected the mergers of ten black
hole binaries via the emission of gravitational radiation. The hypothesis that
these black holes have formed during the cosmic QCD epoch and make up all of
the cosmic dark matter, has been rejected by many authors reasoning that, among
other constraints, primordial black hole (PBH) dark matter would lead to orders
of magnitude larger merger rates than observed. We revisit the calculation of
the present PBH merger rate. Solar mass PBHs form clusters at fairly high
redshifts, which evaporate at lower redshifts. We consider in detail the
evolution of binary properties in such clusters due to three-body interactions
between the two PBH binary members and a third by-passing PBH, for the first
time, by full numerical integration. A Monte-Carlo analysis shows that formerly
predicted merger rates are reduced by orders of magnitude due to such
interactions. The natural prediction of PBH dark matter formed during the QCD
epoch yields a pronounced peak around $1M_{\odot}$ with a small mass fraction
of PBHs on a shoulder around $30M_{\odot}$, dictated by the well-determined
equation of state during the QCD epoch. We employ this fact to make a tentative
prediction of the merger rate of $\sim 30M_{\odot}$ PBH binaries, and find it
very close to that determined by LIGO/Virgo. Furthermore we show that current
LIGO/Virgo limits on the existence of $\sim M_{\odot}$ binaries do not exclude
QCD PBHs to make up all of the cosmic dark matter. Neither do constraints on
QCD PBHs from microlensing, the stochastic gravitational background,
pre-recombination accretion, or dwarf galaxies pose a problem. We caution,
however, in this numerically challenging problem some possibly relevant effects
could not be treated, and merger rates are still somewhat uncertain.