  Biomolecular electrostatics is key in protein function and the chemical
processes affecting it. Implicit-solvent models via the Poisson-Boltzmann (PB)
equation provide insights with less computational cost than atomistic models,
making large-system studies -- at the scale of viruses -- accessible to more
researchers. Here we present a high-productivity and high-performance linear PB
solver based on Exafmm, a fast multipole method library, and Bempp, a Galerkin
boundary element method package. The workflow integrates an easy-to-use Python
interface with optimized computational kernels, and can be run interactively
via Jupyter notebooks, for faster prototyping. Our results show the capability
of the software, confirm code correctness, and assess performance with between
8,000 and 2 million elements. Showcasing the power of this interactive
computing platform, we study the conditioning of two variants of the boundary
integral formulation with just a few lines of code. Mesh-refinement studies
confirm convergence as $1/N$, for $N$ boundary elements, and a comparison with
results from the trusted APBS code using various proteins shows agreement. Our
binding energy calculations using 9 various complexes align with the results
from using five other grid-based PB solvers. Performance results include
timings, breakdowns, and computational complexity. Exafmm offers evaluation
speeds of just a few seconds for tens of millions of points, and
$\mathcal{O}(N)$ scaling. The trend observed in our performance comparison with
APBS demonstrates the advantage of Bempp-Exafmm in applications involving
larger structures or requiring higher accuracy. Computing the solvation free
energy of a Zika virus, represented by 1.6 million atoms and 10 million
boundary elements, took 80-min runtime on a single compute node (dual 20-core).
