Data recorded by the D0 experiment at the Fermilab Tevatron Collider are analyzed to search for neutral Higgs bosons produced in association with $b$ quarks. This production mode can be enhanced in the minimal supersymmetric standard model (MSSM). The search is performed in the three $b$ quark channel using multijet triggered events corresponding to an integrated luminosity of $1 \mathrm{fb}^{-1}$. No statistically significant excess of events with respect to the predicted background is observed and limits are set in the MSSM parameter space.