  The Poisson--Boltzmann equation is widely used to model electrostatics in
molecular systems. Available software packages solve it using finite
difference, finite element, and boundary element methods, where the latter is
attractive due to the accurate representation of the molecular surface and
partial charges, and exact enforcement of the boundary conditions at infinity.
However, the boundary element method is limited to linear equations and
piecewise constant variations of the material properties. In this work, we
present a scheme that couples finite and boundary elements for the
Poisson--Boltzmann equation, where the finite element method is applied in a
confined {\it solute} region, and the boundary element method in the external
{\it solvent} region. As a proof-of-concept exercise, we use the simplest
methods available: Johnson--N\'ed\'elec coupling with mass matrix and diagonal
preconditioning, implemented using the Bempp-cl and FEniCSx libraries via their
Python interfaces. We showcase our implementation by computing the polar
component of the solvation free energy of a set of molecules using a constant
and a Gaussian-varying permittivity. We validate our implementation against the
finite difference code APBS (to 0.5\%), and show scaling from protein G B1 (955
atoms) up to immunoglobulin G (20\,148 atoms). For small problems, the coupled
method was efficient, outperforming a purely boundary integral approach. For
Gaussian-varying permittivities, which are beyond the applicability of boundary
elements alone, we were able to run medium to large sized problems on a single
workstation. Development of better preconditioning techniques and the use of
distributed memory parallelism for larger systems remains an area for future
work. We hope this work will serve as inspiration for future developments for
molecular electrostatics with implicit solvent models.
