Erasing entire clump in simulation

Hi,
I wonder if I am doing correct about deleting clumps in my simulation. My intention is to let the clump fall onto the sphere particles inside the cylinder container but I cut off the code about particles first (it's more than 260K particles) to make sure that the clump would be erase after unbalancedForce is smaller than 0.5. I tried several ways to write it but nothing work, Please give me some piece of advice
Cheers.

####### My code #######
import random
import math
from yade import geom, pack, utils, plot, ymport, export
import numpy as np

# Define cylinder with funnel parameters
center = (0, 0, 0.1)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# add cylinder to simulation
O.bodies.append(cylinder)

# plate properties
clump_plate = CohFrictMat(density=7500, label='clump_plate')

# add properties
O.materials.append(clump_plate)

# clumping parameters
bodyList = []
zi = 0.115
radius_clump = 0.00221
clumpLimit = (diameter/2) - 0.005
# layer and offset condition
num_layers = 8
z_offset_increment = 0.00128

# adjust the z-coordinate based on the layer
for layer in range(num_layers):
    z_offset = layer * z_offset_increment
    # if layer is divisible by 2, move x and y axis by 2mm
    if layer % 2 == 0:
        x_offset = 0.002
        y_offset = 0.002
    else:
        x_offset = 0
        y_offset = 0

    for xi in range(-60, 60, 4):
        x = xi / 1000 + x_offset
        for yi in range(-60, 60, 4):
            y = yi / 1000 + y_offset

            # Check if the sphere is within the clump limit
            if (math.sqrt(x**2 + y**2) + radius_clump) <= clumpLimit:
                sphere_obj = sphere([x, y, zi + z_offset], radius_clump)
                # Assign material properties to the sphere
                sphere_obj.material = clump_plate
                bodyList.append(O.bodies.append(sphere_obj))

# set color
for index, layer in enumerate(O.bodies):
    if not isinstance(layer.shape, Sphere):
        continue
    if index % 2 == 0:
        layer.shape.color = (1, 1, 1) # white color
    else:
        layer.shape.color = (0.5, 0, 0) # navy-blue color

# add clump to bodyList
idClump = O.bodies.clump(bodyList)

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        # handle sphere+sphere and facet+sphere collisions
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
    # call the checkUnbalanced function (defined below) every 600 seconds
    PyRunner(command='checkUnbalanced()', realPeriod=2)
]
O.dt = .5 * PWaveTimeStep()

def checkUnbalanced():

    if unbalancedForce() < .3 :
        # Remove clumps from the simulation
        O.bodies.erase(idClump)
        # Stop the simulation
        O.pause()