I also use fix addforce to apply a force which depends on the z position of the
particles. I use fix ave/spatial to get a number density profile in
the z direction--i.e. rho(z). You already know fext(z), the position
dependent external force, as you supplied it to fix addforce, so you
can get the total external force on all particles by integrating
fext(z)*rho(z) and multiplying by cross-sectional area. All this
assumes that your force is constant with respect to time.