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édélec 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.

What problem does this paper attempt to address?

### What problem does this paper attempt to solve? This paper aims to solve the limitations encountered when using the boundary element method (BEM) to solve the Poisson - Boltzmann equation (PBE) in the electrostatics calculation of molecular solvation. Specifically: 1. **Limitations of the boundary element method**: - BEM is suitable for linear equations and piecewise - constant changes in material properties, but has difficulties in dealing with spatially - varying dielectric constants or nonlinear PBEs. - BEM can only accurately represent molecular surfaces and partial charges, and accurately impose boundary conditions at infinity, but is powerless in the face of complex changes in material properties. 2. **Combining the finite element method (FEM) with the boundary element method (BEM)**: - The paper proposes a method of coupling FEM and BEM to overcome the above - mentioned limitations of BEM. Through this method, FEM is applied to the restricted solute region, while BEM is applied to the external solvent region. - This coupling method can more flexibly handle spatially - varying field parameters, such as Gaussian - varying dielectric constants, which cannot be achieved by pure BEM. 3. **Verification and performance evaluation**: - The paper verifies the effectiveness of its method by calculating the polar solvation free energy of a series of molecules, including cases with constant and Gaussian - varying dielectric constants. - The study also shows the performance of this method on molecules of different scales, from small molecules to large proteins (such as immunoglobulin G), and successfully runs medium - to - large - scale problems on a single workstation. ### Formula summary - **Poisson - Boltzmann equation**: \[ \begin{aligned} &\nabla^2 \phi_1(x)=\frac{1}{\epsilon_1}\sum_{k = 1}^{N_q}q_k\delta(x,x_k),\quad x\in\Omega_1,\\ &(\nabla^2-\kappa^2)\phi_2(x) = 0,\quad x\in\Omega_2,\\ &\phi_1(x)=\phi_2(x),\quad x\in\Gamma,\\ &\epsilon_1\frac{\partial\phi_1}{\partial n}(x)=\epsilon_2\frac{\partial\phi_2}{\partial n}(x),\quad x\in\Gamma. \end{aligned} \] where $\epsilon_1$ and $\epsilon_2$ are the dielectric constants in the solute and solvent respectively, $\kappa$ is the reciprocal of the Debye length, and $q_k$ is the partial charge located at $x_k$. - **Solvation free energy**: \[ \Delta G_{\text{solv}}=\frac{1}{2}\int_{\Omega_1}\rho(x)\phi_r(x)=\frac{1}{2}\sum_{k = 1}^{N_q}q_k\phi_r(x_k). \] - **Gaussian - varying dielectric constant**: \[ \epsilon(r)=(1 - \rho)\epsilon_1+\rho\epsilon_2, \] where \[ \rho(r)=\prod_i\left(1-\exp\left(\frac{\|r - x_i\|}{\sigma^2R_i^2}\right)\right), \] $R_i$ is the van der Waals radius of atom $i$, and $\sigma$ is a scaling factor. ### Conclusion By coupling FEM and BEM, this paper provides an efficient, accurate and flexible computational tool for dealing with molecular electrostatics problems, especially when it comes to spatially - varying field parameters.

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

Accurate boundary-integral formulations for the calculation of electrostatic forces with an implicit-solvent model

Stable non-symmetric coupling of the finite volume and the boundary element method for convection-dominated parabolic-elliptic interface problems

An adaptive non-symmetric finite volume and boundary element coupling method for a fluid mechanics interface problem

Calculation of the Maxwell Stress Tensor and the Poisson-Boltzmann Force on A Solvated Molecular Surface Using Hypersingular Boundary Integrals

Incorporating Dipolar Solvents with Variable Density in Poisson-Boltzmann Electrostatics

The Poisson-Boltzmann model for implicit solvation of electrolyte solutions: Quantum chemical implementation and assessment via Sechenov coefficients

Mixed finite element analysis for a modified Poisson–Fermi interface problem accounting for electrostatic correlations

Domain Decomposition Method for Poisson--Boltzmann Equations based on Solvent Excluded Surface

Accurate, robust and reliable calculations of Poisson-Boltzmann solvation energies

Poisson-Boltzmann model for protein-surface electrostatic interactions and grid-convergence study using the PyGBe code

Optimized Parallelization of Boundary Integral Poisson-Boltzmann Solvers

Some challenges of diffused interfaces in implicit-solvent models

A Nonlinear Boundary Condition for Continuum Models of Biomolecular Electrostatics

An Efficient Finite Element Iterative Method for Solving a Nonuniform Size Modified Poisson-Boltzmann Ion Channel Model

Efficient coupled lattice Boltzmann and Discrete Element Method simulations of small particles in complex geometries

Testing a heterogeneous polarizable continuum model against exact Poisson boundary conditions

A Direct Simulation Approach for the Poisson-Boltzmann Equation Using the Random Batch Method.

Immersed boundary method for dynamic simulation of polarizable colloids of arbitrary shape in explicit ion electrolytes

A Finite Element Solution of Lateral Periodic Poisson–Boltzmann Model for Membrane Channel Proteins