def METHOD_NAME(self):
    RE = self.RE
    N0 = self.N0
    types = self.types
    system = self.system
    # chemical warmup - get close to chemical equilibrium before we start
    # sampling
    RE.reaction(steps=40 * N0)
    average_NH = 0.0
    average_NHA = 0.0
    average_NA = 0.0
    num_samples = 1000
    for _ in range(num_samples):
        RE.reaction(steps=10)
        average_NH += system.number_of_particles(type=types["H+"])
        average_NHA += system.number_of_particles(type=types["HA"])
        average_NA += system.number_of_particles(type=types["A-"])
    average_NH /= float(num_samples)
    average_NA /= float(num_samples)
    average_NHA /= float(num_samples)
    average_alpha = average_NA / float(N0)
    # note you cannot calculate the pH via -log10(<NH>/volume) in the
    # constant pH ensemble, since the volume is totally arbitrary and does
    # not influence the average number of protons
    pH = self.pH
    pKa = self.pKa
    target_alpha = Test.ideal_alpha(pH)
    rel_error_alpha = abs(average_alpha - target_alpha) / target_alpha
    # relative error
    self.assertLess(
        rel_error_alpha,
        0.015,
        msg="\nDeviation from ideal titration curve is too big for the given input parameters.\n"
        + f"  pH: {pH:.2f}"
        + f"  pKa: {pKa:.2f}"
        + f"  average NH: {average_NH:.1f}"
        + f"  average NA: {average_NA:.1f}"
        + f"  average NHA: {average_NHA:.1f}"
        + f"  average alpha: {average_alpha:.3f}"
        + f"  target alpha: {target_alpha:.3f}"
        + f"  rel_error: {rel_error_alpha * 100:.1f}%"
    )