How to get the resultant force on a particle

Hi,

I'd like to ask that how to get the resultant force on a particle.

I came out with two ways to get the (magnitude of) resultant force on a particle as follow.

(1) using O.forces.f(id)
(2) calculating from the interactions on that particle:
a=[]
for i in O.bodies[id].intrs():
        f = i.phys.normalForce + i.phys.shearForce
        a.append(f)

resultF=Vector3(0,0,0)
for i in range(len(a)):
        resultF=resultF+a[i]

However, in some of situations, these two ways can have the same magnitude of resultant force, for instance, using MWE1, the magnitude of resultant force obtained from the two approaches are the same. But using MWE2, I got different results from the two approaches. (both the two MWE are slightly modified from YADE example[1])

Could you please help me with: whether the way I calculate the resultant force is correct; if they are correct, why we can get different results?
########### MWE1 which can lead to the same results###########
O.bodies.append(
        [
                # fixed: particle's position in space will not change (support)
                sphere(center=(0, 0, 0), radius=.5, fixed=True),
                # this particles is free, subject to dynamics
                sphere((0, 0, 0.8), .5),
                sphere((0, 0.8, 0), .5)
        ]
)

# FUNCTIONAL COMPONENTS

# simulation loop -- see presentation for the explanation
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom()], # collision geometry
                [Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
                [Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
        ),
        # Apply gravity force to particles. damping: numerical dissipation of energy.
        NewtonIntegrator(gravity=(0, 0, 0), damping=0.1)
]

# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt = .5e-4 * PWaveTimeStep()

O.step()

print O.forces.f(0)

a=[]
for i in O.bodies[0].intrs():
        f = i.phys.normalForce + i.phys.shearForce
        a.append(f)

resultF=Vector3(0,0,0)
for i in range(len(a)):
        resultF=resultF+a[i]

print resultF

########### MWE2 which can lead to different results###########

readParamsFromTable(rMean=.05, rRelFuzz=.3, maxLoad=1e6, minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=rMean, rRelFuzz=rRelFuzz)
sp.toSimulation()

O.engines = [
        ForceResetter(),
        # sphere, facet, wall
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(
                # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
        # the label creates an automatic variable referring to this engine
        # we use it below to change its attributes from the functions called
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = .5 * PWaveTimeStep()

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
        # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
        if O.iter < 5000:
                return
        # the rest will be run only if unbalanced is < .1 (stabilized packing)
        if unbalancedForce() > .1:
                return
        # add plate at the position on the top of the packing
        # the maximum finds the z-coordinate of the top of the topmost particle
        O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
        global plate # without this line, the plate variable would only exist inside this function
        plate = O.bodies[-1] # the last particles is the plate
        # Wall objects are "fixed" by default, i.e. not subject to forces
        # prescribing a velocity will therefore make it move at constant velocity (downwards)
        plate.state.vel = (0, 0, -.1)
        # start plotting the data now, it was not interesting before
        O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
        # next time, do not call this function anymore, but the next one (unloadPlate) instead
        checker.command = 'unloadPlate()'

def unloadPlate():
        # if the force on plate exceeds maximum load, start unloading
        if abs(O.forces.f(plate.id)[2]) > maxLoad:
                plate.state.vel *= -1
                # next time, do not call this function anymore, but the next one (stopUnloading) instead
                checker.command = 'stopUnloading()'

def stopUnloading():
        if abs(O.forces.f(plate.id)[2]) < minLoad:
                # O.tags can be used to retrieve unique identifiers of the simulation
                # if running in batch, subsequent simulation would overwrite each other's output files otherwise
                # d (or description) is simulation description (composed of parameter values)
                # while the id is composed of time and process number
                plot.saveDataTxt(O.tags['d.id'] + '.txt')
                O.pause()

def addPlotData():
        if not isinstance(O.bodies[-1].shape, Wall):
                plot.addData()
                return
        Fz = O.forces.f(plate.id)[2]
        plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)

O.run(10000,1)

waitIfBatch()

print O.forces.f(20)

a=[]
for i in O.bodies[20].intrs():
        f = i.phys.normalForce + i.phys.shearForce
        a.append(f)

resultF=Vector3(0,0,0)
for i in range(len(a)):
        resultF=resultF+a[i]

print resultF

Thanks
Leonard

[1]https://yade-dev.gitlab.io/trunk/tutorial-examples.html