def METHOD_NAME(self, T, reaction, forward=True):
    """
    Return the diffusive limit on the rate coefficient, k_diff.
    This is the upper limit on the rate, in the specified direction.
    (ie. forward direction if forward=True [default] or reverse if forward=False)
    Returns the rate coefficient k_diff in m3/mol/s.
    """
    if forward:
        reacting = reaction.reactants
    else:
        reacting = reaction.products
    if len(reacting) < 2:
        raise Exception("Cannot calculate diffusion limit for a unimolecular reaction")
    radii = 0.0
    diffusivities = []
    for spec in reacting:
        solute_data = self.database.get_solute_data(spec)
        # calculate radius with the McGowan volume and assuming sphere
        radius = ((75 * solute_data.V / constants.pi / constants.Na) ** (1. / 3)) / 100  # m
        diff = solute_data.get_stokes_diffusivity(T, self.get_solvent_viscosity(T))
        radii += radius  # m
        diffusivities.append(diff)  # m^2/s
    # Calculate Smoluchowski kinetics for any reaction order
    # Flegg, SIAM J. Appl. Math., 76 (4), 2016
    N = len(reacting)
    Dinv = 1.0 / np.array(diffusivities)
    Dhat = np.empty(N - 1)
    alpha = (3.0 * N - 5.0) / 2.0
    denom = 0.0
    for i in range(N - 1):
        Dbar = 1.0 / np.sum(Dinv[:(i + 1)])
        Dhat[i] = diffusivities[i + 1] + Dbar
        for j in range(i + 1, N):
            denom += 1.0 / (diffusivities[i] * diffusivities[j])
    delta = np.sum(Dinv) / denom
    k_diff = (np.prod(Dhat) ** 1.5
              * 4.0 * constants.pi ** (alpha + 1.0) / math.gamma(alpha)
              * (radii / np.sqrt(delta)) ** (2.0 * alpha)
              * constants.Na ** (N - 1.0))  # m^(3*(N-1))/mol^(N-1)/s
    return k_diff