Just to clarify:  are you using APBS for the solvation energy and some
other program ( e.g., the provided `coulomb` program) for the Coulombic
energies?

I use APBS for the solvation energy and the 'coulomb.py' script for the Coulombic energies.

> Then would it be accurate to say that the effective force at r = (x,y,z) is
> (delta E)/(delta r), where E is the solvation energy plus coulombic energy?

Are you referring to the analytic forces returned by APBS or the
forces you get by evaluating energies at slightly different positions
and subtracting the results?

The latter.
 
Thanks,

Bihui