def METHOD_NAME(simulation_factory, lattice_snapshot_factory, phi, eta_p,
                ntrial, use_clusters):
    # number of spheres
    n = 5
    N = n**3
    d_sphere = 1.0
    V_sphere = math.pi / 6.0 * math.pow(d_sphere, 3.0)
    # depletant-colloid size ratio (unity, to emphasize many-body effects)
    q = 1.0
    L_target = math.pow(N * V_sphere / phi, 1.0 / 3.0)
    a = L_target / n
    sim = simulation_factory(
        lattice_snapshot_factory(particle_types=['A', 'B'],
                                 dimensions=3,
                                 a=a,
                                 n=n))
    mc = hoomd.hpmc.integrate.Sphere(seed=1, default_d=0.1, default_a=0.1)
    sim.operations.integrator = mc
    mc.shape['A'] = dict(diameter=d_sphere)
    mc.shape['B'] = dict(diameter=d_sphere * q)
    tune = hoomd.hpmc.tune.MoveSize.scale_solver(trigger=hoomd.trigger.And(
        [hoomd.trigger.Periodic(100),
         hoomd.trigger.Before(1000)]),
                                                 moves=['d'],
                                                 types=['A'],
                                                 max_translation_move=d_sphere,
                                                 gamma=1,
                                                 target=0.2)
    sim.operations.tuners.append(tune)
    sim.run(1000)
    # warm up
    sim.run(2000)
    # set depletant fugacity
    nR = eta_p / (math.pi / 6.0 * math.pow(d_sphere * q, 3.0))
    mc.depletant_fugacity[('B', 'B')] = nR
    mc.depletant_ntrial[('B', 'B')] = ntrial
    if use_clusters:
        cl = hoomd.hpmc.update.Clusters(trigger=1, seed=1)
        sim.operations.updaters.append(cl)
        # use clusters exclusively to equilibrate
        mc.d['A'] = 0
    sim.run(100)