def METHOD_NAME(self):
    # reactant and product types
    type_A = 0
    type_B = 1
    type_C = 2
    type_D = 3
    type_E = 4
    types = [type_A, type_B, type_C, type_D, type_E]
    # reactant and product stoichiometric coefficients
    nu_A = -2
    nu_B = -3
    nu_C = 4
    nu_D = 1
    nu_E = 3
    nu_s = [nu_A, nu_B, nu_C, nu_D, nu_E]
    # reaction constant
    K = 0.001
    # reference concentration: 1 mol/l
    c_ref = 1.0
    # simulation units: 1 sigma = 3.55 Angstrom
    conv_sim_unit_to_mol_per_l = 37.1
    Gamma = K * (c_ref / conv_sim_unit_to_mol_per_l)**np.sum(nu_s)
    volume = self.system.volume()
    # setup N0 batches of reactants
    N0 = 30
    self.system.part.add(
        pos=np.random.random((-nu_A * N0, 3)) * self.system.box_l,
        type=-nu_A * N0 * [type_A])
    self.system.part.add(
        pos=np.random.random((-nu_B * N0, 3)) * self.system.box_l,
        type=-nu_B * N0 * [type_B])
    # use an exclusion radius of 0 to simulate an ideal gas
    RE = espressomd.reaction_methods.ReactionEnsemble(
        kT=1., exclusion_range=0., seed=42)
    RE.add_reaction(
        gamma=Gamma, reactant_types=[type_A, type_B],
        reactant_coefficients=[abs(nu_A), abs(nu_B)],
        product_types=[type_C, type_D, type_E],
        product_coefficients=[nu_C, nu_D, nu_E],
        default_charges={type_A: 0, type_B: 0, type_C: 0, type_D: 0, type_E: 0})
    # Set the hidden particle type to the lowest possible number to speed
    # up the simulation
    RE.set_non_interacting_type(type=max(types) + 1)
    # warmup
    RE.reaction(steps=50)
    numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []}
    for _ in range(40):
        RE.reaction(steps=6)
        for key in types:
            numbers[key].append(self.system.number_of_particles(type=key))
    def calculate_K(c_s, nu_s):
        return np.prod(
            [(c_i / c_ref)**nu_i for c_i, nu_i in zip(c_s, nu_s)])
    def equations(concentrations):
        c_A, c_B, c_C, c_D, c_E = concentrations
        eq1 = K - calculate_K(concentrations, nu_s)
        eq2 = N0 - (1.0 / abs(nu_A) * c_A / conv_sim_unit_to_mol_per_l +
                    1.0 / nu_D * c_D / conv_sim_unit_to_mol_per_l) * volume
        eq3 = c_A / c_B - float(nu_A) / float(nu_B)
        eq4 = c_C / c_D - float(nu_C) / float(nu_D)
        eq5 = c_C / c_E - float(nu_C) / float(nu_E)
        return (eq1, eq2, eq3, eq4, eq5)
    concentrations_sim = {
        type_i: np.mean(n_i) / volume * conv_sim_unit_to_mol_per_l
        for type_i, n_i in numbers.items()}
    concentrations_analytic = dict(zip(types, scipy.optimize.fsolve(
        equations, list(concentrations_sim.values()))))
    K_sim = calculate_K(concentrations_sim.values(), nu_s)
    N0_sim = (1.0 / abs(nu_A) * concentrations_sim[type_A] + 1.0 / nu_D *
              concentrations_sim[type_D]) / conv_sim_unit_to_mol_per_l * volume
    err_msg = "Concentration of species {} doesn't match analytical result"
    for key in types:
        self.assertAlmostEqual(concentrations_sim[key],
                               concentrations_analytic[key],
                               delta=1e-3,
                               msg=err_msg.format(key))
    self.assertAlmostEqual(K_sim, K, delta=1e-3)
    self.assertAlmostEqual(N0_sim, N0, delta=1e-3)