I also tried to run it with only 100 and 1000 events and in these cases it works without any problem, but with 10000 I get the same error than previously. I have also checked that I have enough disk memory to run these simulations (especially on the clusters where hundreds of GB are available) so the problem shouldn't be that the files takes too much space when I increase the number of events.

With that line, the computation is going trough but I see some Fortran warning:
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_OVERFLOW_FLAG
Which seems to indicate that the model has some division by zero.

I will investigate those and report those to the author of the model if needed.

Actually, the numerical issue is due to an inconsistent setup on your side.
You are using a 4 flavour model to do a 5 flavour computation.

In addition to not be advised, you do not set correctly MG5aMC to do a 5 flavour computation since you simple use default run_card
which is set for a 4 flavour and not for a 5 flavour one (just setting the b mass to zero is not enough to have the correct result).

This lead to a problem with the scale computation which is not correctly defined and raises series problems including the one with the discrete sampler pointed above.

Thank you for your answer! I actually wanted to do the calculation in the 4 flavour scheme but I forgot to set the mass of the b quark to its usual value. I changed that and everything is now running without any problem.