def METHOD_NAME(self, single_precision: bool):
    relative_precision: float = 1E-6 if single_precision else 1E-7
    lattice = espressomd.electrokinetics.LatticeWalberla(
        n_ghost_layers=1, agrid=self.AGRID)
    eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)
    self.system.ekcontainer = espressomd.electrokinetics.EKContainer(
        tau=self.TAU, solver=eksolver)
    reaction_rate: float = 1e-5
    stoech_coeffs = [2.0, 1.0, 1.2, 2.2]
    product_coeff = 1.5
    educt_species = []
    reactants = []
    for coeff in stoech_coeffs:
        species = espressomd.electrokinetics.EKSpecies(
            lattice=lattice, density=coeff * self.INITIAL_DENSITY,
            diffusion=self.DIFFUSION_COEFFICIENT, valency=0.0,
            advection=False, friction_coupling=False,
            single_precision=single_precision, tau=self.TAU)
        self.system.ekcontainer.add(species)
        reactants.append(
            espressomd.electrokinetics.EKReactant(
                ekspecies=species,
                stoech_coeff=-coeff,
                order=coeff))
        educt_species.append(species)
    ek_species_product = espressomd.electrokinetics.EKSpecies(
        lattice=lattice, density=0.0, diffusion=self.DIFFUSION_COEFFICIENT,
        valency=0.0, advection=False, friction_coupling=False,
        single_precision=single_precision, tau=self.TAU)
    self.system.ekcontainer.add(ek_species_product)
    reactants.append(
        espressomd.electrokinetics.EKReactant(
            ekspecies=ek_species_product,
            stoech_coeff=product_coeff,
            order=0.0))
    reaction = espressomd.electrokinetics.EKBulkReaction(
        reactants=reactants, coefficient=reaction_rate, lattice=lattice, tau=self.TAU)
    self.system.ekcontainer.reactions.add(reaction)
    self.assertEqual(len(self.system.ekcontainer.reactions), 1)
    self.system.integrator.run(self.TIME)
    domain_volume = np.product(ek_species_product.shape)
    analytic_time = (self.TIME + 0.5) * self.system.time_step
    measured_educt_densities = np.zeros(len(stoech_coeffs))
    for i, educt in enumerate(educt_species):
        measured_educt_densities[i] = np.sum(
            educt[:, :, :].density) / domain_volume
    measured_product_density = np.sum(
        ek_species_product[:, :, :].density) / domain_volume
    analytic_educt_densities = np.zeros(len(stoech_coeffs))
    for i, coeff in enumerate(stoech_coeffs):
        analytic_educt_densities[i] = coeff * self.analytic_density_base(
            analytic_time, stoech_coeffs, reaction_rate, self.INITIAL_DENSITY)
    analytic_product_density = product_coeff * \
        (self.INITIAL_DENSITY -
            self.analytic_density_base(analytic_time, stoech_coeffs,
                                       reaction_rate, self.INITIAL_DENSITY))
    np.testing.assert_allclose(
        measured_educt_densities / measured_educt_densities[0],
        analytic_educt_densities / analytic_educt_densities[0],
        rtol=relative_precision, atol=0)
    np.testing.assert_allclose(
        measured_educt_densities,
        analytic_educt_densities,
        rtol=2 * reaction_rate,
        atol=0)
    np.testing.assert_allclose(
        measured_product_density,
        analytic_product_density,
        rtol=2 *
        reaction_rate *
        len(stoech_coeffs),
        atol=0)
    self.system.ekcontainer.reactions.remove(reaction)
    self.assertEqual(len(self.system.ekcontainer.reactions), 0)