Coupling finite and boundary element methods to solve the Poisson--Boltzmann equation for electrostatics in molecular solvation

TL;DR

This paper introduces a coupled finite and boundary element method for solving the Poisson--Boltzmann equation, enabling accurate electrostatic modeling in molecular solvation with improved efficiency and applicability to complex permittivity variations.

Contribution

It presents a novel coupling scheme combining finite and boundary element methods for the Poisson--Boltzmann equation, extending applicability to Gaussian-varying permittivities and larger molecular systems.

Findings

Validated against APBS with 0.5% accuracy

Outperformed pure boundary element methods on small problems

Enabled medium to large problem solving on a single workstation

Abstract

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.