The grand canonical ensemble describes a system at fixed T and chemical potential μ exchanging particles with a reservoir. The mean occupation follows the Fermi–Dirac or Bose–Einstein distribution. Here we simulate a multilevel system at chemical equilibrium.